Logic

Salvage topsoil was removed from a donor site and delivered to three recipient sites in 2015. We sampled all sites (Donor and recipients) prior to topsoil delivery and sequenced AM Fungi from all sites in 2015 and 2017 from manipulative plots, set up in a randomized block design at each site with: control treatments (no topsoil added), a dusting of topsoil, and three levels of topsoil thicknesses (2“, 4” and 6" thick layers of this topsoil, originating from the donor site)), delivered to all recipient sites).

knitr::opts_chunk$set(cache=TRUE, tidy=TRUE, error = FALSE, eval = TRUE, message = FALSE, warning = FALSE, rows.print=5, cols.min.print=4, fig.width=6, fig.height=4.5)

Questions:

**Q1: Propagule determination: Do AM fungal communities resemble the donor site, regardless of where soil was delivered? Is provenance a driver of AM fungal community composition?

**Q2: Environmental filtering: Do AM fungal communities resemble the recipient sites, and are more dissimilar to the AM fungal communities from the donor site?

**Q3: Propagule pressure: Do AM fungal communities resemble the donor site more in higher level topsoil treatments, than they do in the control or dusted sites?

Outputs:

Analysis of similarity:

Determine the similarity of AMF communities at each recipient site to the donor site, from where topsoils originated from. Generate heatmaps with distances, to visually determine - which sites/samples are more similar to each other, and which are more dissimilar. Generate NMDS / Pcoa, plotted with Year as shapes and Sites as colors, and Description as …? Run permanova Do AM fungal communities at each site resemble the recipient sites, and are they more dissimilar to the donor sites?

Functional group richness:

Determine the richness of AMF functional groups in the donor site, from where topsoils originated from, and in the recipient sites pre-topsoil delivery (pre-treatment), and post-topsoil delivery (post-treatment)

Taxonomic diversity:

Determine the OTU taxa richness (alpha diversity) of AMF communities in the donor site, from where topsoils originated from, and in the recipient sites pre-topsoil delivery (pre-treatment), and post-topsoil delivery (post-treatment). Determine the beta diversity of AMF communities in the donor site, from where topsoils originated from, and in the recipient sites pre-topsoil delivery (pre-treatment), and post-topsoil delivery (post-treatment).

Topsoil level treatments

Propagule pressure: Do AM fungal communities resemble the donor site more in thicker topsoil layer treatments than they do in the control or dusted sites? Are the sites with Treat groups = S, more similar to the donor site than F, and are treatment groups S and F more similar to the donor site than T? Is this relationship clinal, with S being the thickest and most similar to donor site, followed by F, and then by T. TO FIGURE OUT - how to address this question statistically??multiple regression with Treat on the X and similarity to Donor Site on the Y?? A permanova with multiple comparisons?

knitr::opts_chunk$set(cache = TRUE, tidy = TRUE, error = FALSE, eval = TRUE, 
    message = FALSE, warning = FALSE, rows.print = 5, cols.min.print = 4, fig.width = 6, 
    fig.height = 4.5)

Load Required Packages

library(tidyverse)
library(dplyr)  ## for data wrangling - %>% function
library(reshape2)  ##melt and cast data
library(ggplot2)  # plotting
library(data.table)
library(stringr)
library(tidyr)  # 'separate' function
library(readxl)  #read xlsx files into r on mac computer
library(vegan)  # dissimilarity matrix, permanova functions
library(magrittr)
library(cowplot)
library(formatR)


date <- format(Sys.time(), "%Y%b%d")
date
## [1] "2018Sep28"

** the following is adapted from 2017-07-07_Chapter_2_SSU_IntialClean **

Import data files

# read in OTU table
ssu <- fread("CSS_table_sorted.txt")
str(ssu)
## Classes 'data.table' and 'data.frame':   386 obs. of  63 variables:
##  $ #OTU ID : chr  "denovo538" "denovo1813859" "denovo1031897" "denovo34" ...
##  $ MM933M  : num  5.29 4.33 7.32 0 0 ...
##  $ MM936M  : num  6.27 11.24 3.4 0 0 ...
##  $ MM941M  : num  4.4 2.95 0 0 0 ...
##  $ MM943M  : num  5.31 4.04 0 0 0 ...
##  $ MM9362CU: num  3.13 3.13 8.77 0 1.56 ...
##  $ MMS700  : num  2.5 6.82 7.82 0 2.5 ...
##  $ MMS701  : num  1.87 2.66 8.27 2.66 0 ...
##  $ MMS702  : num  0 2.79 8.24 2.79 0 ...
##  $ MMS703  : num  0 2.58 7.18 3.45 0 ...
##  $ MMS704  : num  2.75 0 7.43 9.75 0 ...
##  $ MMS705  : num  3.28 0 6.99 0 4.76 ...
##  $ MMS706  : num  1.83 2.61 8.32 0 0 ...
##  $ MMS707  : num  3.03 7.49 7.43 12.04 0 ...
##  $ MMS708  : num  0 3.88 7.64 2.97 0 ...
##  $ MMS709  : num  5.69 8.79 7.95 0 1.56 ...
##  $ MMS710  : num  3.98 3.07 8.5 2.23 0 ...
##  $ MMS711  : num  4.44 1.99 8.14 1.99 0 ...
##  $ MMS712  : num  2.23 3.07 5.91 7.39 4.75 ...
##  $ MMS713  : num  2.87 2.05 8.23 2.05 0 ...
##  $ MMS714  : num  4.03 4.34 6.54 2.28 8.78 ...
##  $ MMS715  : num  3.03 6.25 4.71 9.76 3.03 ...
##  $ MMS716  : num  4.22 2.43 8.93 10.72 0 ...
##  $ MMS717  : num  3.12 5.56 8.23 2.28 0 ...
##  $ MMS890  : num  4.73 2.89 9.24 0 0 ...
##  $ MMS891  : num  4.35 2.55 8.66 0 0 ...
##  $ MMS900  : num  0 2.95 7.12 2.95 0 ...
##  $ MMS901  : num  2.56 6.47 6.11 2.56 0 ...
##  $ MMS902  : num  4.22 7.82 6.94 7.64 0 ...
##  $ MMS903  : num  4.52 6.31 9.14 2.86 0 ...
##  $ MMS904  : num  3.7 2.32 9.18 2.32 0 ...
##  $ MMS905  : num  2.95 6.57 7.42 0 0 ...
##  $ MMS906  : num  3.5 3.5 9.15 0 0 ...
##  $ MMS907  : num  3.54 4.88 7.56 4.69 0 ...
##  $ MMS908  : num  5.26 6.77 9 3.22 0 ...
##  $ MMS909  : num  3.76 0 6.25 0 0 ...
##  $ MMS910  : num  3.4 0 7.02 6.86 0 ...
##  $ MMS911  : num  0 0 7.01 4.39 5.73 ...
##  $ MMS912  : num  3.01 2.18 8.48 0 0 ...
##  $ MMS913  : num  3.94 1.95 9.29 1.95 0 ...
##  $ MMS914  : num  3.83 0 7.66 2.44 4.53 ...
##  $ MMS915  : num  0 3.72 6.21 2.02 0 ...
##  $ MMS916  : num  4.04 2.03 7.68 4.69 0 ...
##  $ MMS917  : num  2.58 0 0 0 7.61 ...
##  $ MMS918  : num  2.78 8.18 0 0 0 ...
##  $ MMS919  : num  3.79 2.4 8.47 0 0 ...
##  $ MMS920  : num  4.31 2.26 6.97 3.1 0 ...
##  $ MMS921  : num  4.16 2.6 7.69 3.68 1.43 ...
##  $ MMS922  : num  3.6 1.92 7.99 2.71 1.92 ...
##  $ MMS923  : num  3.55 3.03 8.85 3.55 0 ...
##  $ MMS924  : num  3.68 0 8.09 2.79 0 ...
##  $ MMS925  : num  3.69 2.31 7.78 3.16 0 ...
##  $ MMS928  : num  3.4 3.4 7.87 2.53 0 ...
##  $ MMS929  : num  4.19 4.89 9.17 1.95 0 ...
##  $ MMS930  : num  4.75 3.81 9.89 0 0 ...
##  $ MMS931  : num  0 4.61 8.61 0 0 ...
##  $ MMS932  : num  4.63 0 9.28 0 0 ...
##  $ MMS938  : num  2.5 2.5 8.11 0 0 ...
##  $ MMS939  : num  6.1 3.73 8.08 0 0 ...
##  $ MMS942  : num  0 0 8.1 2.97 0 ...
##  $ MMS943  : num  0 2.91 7.07 4.92 0 ...
##  $ MMS944  : num  5.13 3.87 8.23 2.96 0 ...
##  $ taxonomy: chr  "k__Fungi; p__Glomeromycota; c__Glomeromycetes; o__Glomerales; f__Glomeraceae; g__Glomus; s__NES27" "k__Fungi; p__Glomeromycota; c__Glomeromycetes; o__Glomerales; f__Glomeraceae; g__Glomus; s__MH3_B" "No blast hit" "k__Fungi; p__Glomeromycota; c__Glomeromycetes; o__Glomerales; f__Glomeraceae; g__Glomus; s__NES27" ...
##  - attr(*, ".internal.selfref")=<externalptr>
## mapping data
metassu <- fread("mapSal18927.txt", header = TRUE)
metassu$Year <- as.character(metassu$Year)
metassu$Rep <- as.factor(metassu$Rep)
metassu$Treat <- as.factor(metassu$Treat)
metassu$HL <- as.numeric(metassu$HL)
str(metassu)
## Classes 'data.table' and 'data.frame':   61 obs. of  13 variables:
##  $ #SampleID           : int  34 35 36 41 46 42 43 6 37 38 ...
##  $ BarcodeSequence     : chr  "GAGATCGTTGGGTGAT" "GAGATCGTTACTGTGC" "GAGATCGTAGAGTGTG" "ATGCTTGGTGGTAACC" ...
##  $ LinkerPrimerSequence: logi  NA NA NA NA NA NA ...
##  $ Location            : chr  "B5" "C5" "D5" "A6" ...
##  $ Project             : chr  "MM-Salvage" "MM-Salvage" "MM-Salvage" "MM-Salvage" ...
##  $ Year                : chr  "2015" "2015" "2015" "2015" ...
##  $ Site                : chr  "DS" "DS" "DS" "DS" ...
##  $ Treat               : Factor w/ 6 levels "C","D","F","O",..: 4 4 4 4 4 4 4 4 4 4 ...
##  $ Rep                 : Factor w/ 6 levels "1","2","3","4",..: 1 2 3 4 5 1 2 1 2 3 ...
##  $ HL                  : num  1.673 1.519 0.669 0.519 0.679 ...
##  $ TSDel               : chr  NA NA NA NA ...
##  $ Description         : chr  "Donor" "Donor" "Donor" "Donor" ...
##  $ SampleID2           : chr  "MMS938" "MMS939" "MMS890" "MMS891" ...
##  - attr(*, ".internal.selfref")=<externalptr>

Clean SSU headers

colnames(ssu)
##  [1] "#OTU ID"  "MM933M"   "MM936M"   "MM941M"   "MM943M"   "MM9362CU"
##  [7] "MMS700"   "MMS701"   "MMS702"   "MMS703"   "MMS704"   "MMS705"  
## [13] "MMS706"   "MMS707"   "MMS708"   "MMS709"   "MMS710"   "MMS711"  
## [19] "MMS712"   "MMS713"   "MMS714"   "MMS715"   "MMS716"   "MMS717"  
## [25] "MMS890"   "MMS891"   "MMS900"   "MMS901"   "MMS902"   "MMS903"  
## [31] "MMS904"   "MMS905"   "MMS906"   "MMS907"   "MMS908"   "MMS909"  
## [37] "MMS910"   "MMS911"   "MMS912"   "MMS913"   "MMS914"   "MMS915"  
## [43] "MMS916"   "MMS917"   "MMS918"   "MMS919"   "MMS920"   "MMS921"  
## [49] "MMS922"   "MMS923"   "MMS924"   "MMS925"   "MMS928"   "MMS929"  
## [55] "MMS930"   "MMS931"   "MMS932"   "MMS938"   "MMS939"   "MMS942"  
## [61] "MMS943"   "MMS944"   "taxonomy"
head(ssu$taxonomy)
## [1] "k__Fungi; p__Glomeromycota; c__Glomeromycetes; o__Glomerales; f__Glomeraceae; g__Glomus; s__NES27"               
## [2] "k__Fungi; p__Glomeromycota; c__Glomeromycetes; o__Glomerales; f__Glomeraceae; g__Glomus; s__MH3_B"               
## [3] "No blast hit"                                                                                                    
## [4] "k__Fungi; p__Glomeromycota; c__Glomeromycetes; o__Glomerales; f__Glomeraceae; g__Glomus; s__NES27"               
## [5] "k__Fungi; p__Glomeromycota; c__Glomeromycetes; o__Glomerales; f__Glomeraceae; g__Glomus; s__sp"                  
## [6] "k__Fungi; p__Glomeromycota; c__Glomeromycetes; o__Diversisporales; f__Acaulosporaceae; g__Acaulospora; s__Acau16"
# rename columns
names(ssu)[1] <- "ssuotu"  #rename first column
names(ssu)[names(ssu) == "taxonomy"] <- "ssutaxonomy"  # rename column that is currently called taxonomy

# split taxonomy column
`?`(str_match)
ssu$ssukingdom <- str_match(ssu$ssutaxonomy, "k__(.*?);")[, 2]
ssu$ssuphylum <- str_match(ssu$ssutaxonomy, "p__(.*?);")[, 2]
ssu$ssuclass <- str_match(ssu$ssutaxonomy, "c__(.*?);")[, 2]
ssu$ssuorder <- str_match(ssu$ssutaxonomy, "o__(.*?);")[, 2]
ssu$ssufamily <- str_match(ssu$ssutaxonomy, "f__(.*?);")[, 2]
ssu$ssugenus <- str_match(ssu$ssutaxonomy, "g__(.*?);")[, 2]
ssu$ssuspecies <- str_match(ssu$ssutaxonomy, "s__(.*?)$")[, 2]

unique(ssu$ssuspecies)
##  [1] "NES27"                               
##  [2] "MH3_B"                               
##  [3] NA                                    
##  [4] "sp"                                  
##  [5] "Acau16"                              
##  [6] "MO_G47"                              
##  [7] "laccatum"                            
##  [8] "acnaGlo2"                            
##  [9] "MO_Ar1"                              
## [10] "leptoticha"                          
## [11] "Whitfield_type_7"                    
## [12] "Torrecillas12b_Glo_G5"               
## [13] "Sanchez_Castro12b_GLO12"             
## [14] "Winther07_D"                         
## [15] "MO_G18"                              
## [16] "intraradices"                        
## [17] "Alguacil12a_Para_1"                  
## [18] "Douhan9"                             
## [19] "Alguacil12b_GLO_G3"                  
## [20] "Glo39"                               
## [21] "Alguacil10_Glo1"                     
## [22] "Ligrone07_sp"                        
## [23] "PT6"                                 
## [24] "MO_G41"                              
## [25] "MO_GB1"                              
## [26] "A1"                                  
## [27] "Liu2012b_Phylo_5"                    
## [28] "Yamato09_C1"                         
## [29] "acnaGlo7"                            
## [30] "MO_G7"                               
## [31] "Alguacil10_Glo6"                     
## [32] "Alguacil12b_PARA1"                   
## [33] "Glo58"                               
## [34] "lamellosum"                          
## [35] "VeGlo18"                             
## [36] "Liu2012a_Ar_1"                       
## [37] "Div"                                 
## [38] "fennica"                             
## [39] "Voyriella_parviflora_symbiont_type_1"
## [40] "Para1_OTU2"                          
## [41] "Whitfield_type_3"                    
## [42] "Aca"                                 
## [43] "MO_G20"                              
## [44] "pyriformis"                          
## [45] "Glo7"                                
## [46] "MO_GC1"                              
## [47] "Winther07_B"                         
## [48] "mosseae"                             
## [49] "caledonium"                          
## [50] "brasilianum"                         
## [51] "Alguacil09b_Glo_G8"                  
## [52] "Schechter08_Arch1"                   
## [53] "Wirsel_OTU21"                        
## [54] "Glo71"                               
## [55] "spurca"                              
## [56] "Glo59"                               
## [57] "Wirsel_OTU16"                        
## [58] "Glo32"                               
## [59] "MO_A8"                               
## [60] "MO_G27"                              
## [61] "MO_G5"
colnames(ssu)
##  [1] "ssuotu"      "MM933M"      "MM936M"      "MM941M"      "MM943M"     
##  [6] "MM9362CU"    "MMS700"      "MMS701"      "MMS702"      "MMS703"     
## [11] "MMS704"      "MMS705"      "MMS706"      "MMS707"      "MMS708"     
## [16] "MMS709"      "MMS710"      "MMS711"      "MMS712"      "MMS713"     
## [21] "MMS714"      "MMS715"      "MMS716"      "MMS717"      "MMS890"     
## [26] "MMS891"      "MMS900"      "MMS901"      "MMS902"      "MMS903"     
## [31] "MMS904"      "MMS905"      "MMS906"      "MMS907"      "MMS908"     
## [36] "MMS909"      "MMS910"      "MMS911"      "MMS912"      "MMS913"     
## [41] "MMS914"      "MMS915"      "MMS916"      "MMS917"      "MMS918"     
## [46] "MMS919"      "MMS920"      "MMS921"      "MMS922"      "MMS923"     
## [51] "MMS924"      "MMS925"      "MMS928"      "MMS929"      "MMS930"     
## [56] "MMS931"      "MMS932"      "MMS938"      "MMS939"      "MMS942"     
## [61] "MMS943"      "MMS944"      "ssutaxonomy" "ssukingdom"  "ssuphylum"  
## [66] "ssuclass"    "ssuorder"    "ssufamily"   "ssugenus"    "ssuspecies"
ssu.save <- ssu
# View(ssu)

# remove Geosiphonaceae
unique(ssu$ssufamily)
## [1] "Glomeraceae"          NA                     "Acaulosporaceae"     
## [4] "Paraglomeraceae"      "Archaeosporaceae"     "Ambisporaceae"       
## [7] "Claroideoglomeraceae" "Diversisporaceae"     "Geosiphonaceae"
drop.family <- c("Geosiphonaceae")
ssu <- ssu[-which(ssu$ssufamily %in% drop.family), ]



# remove no blast hits
unique(ssu$ssutaxonomy)
##  [1] "k__Fungi; p__Glomeromycota; c__Glomeromycetes; o__Glomerales; f__Glomeraceae; g__Glomus; s__NES27"                                  
##  [2] "k__Fungi; p__Glomeromycota; c__Glomeromycetes; o__Glomerales; f__Glomeraceae; g__Glomus; s__MH3_B"                                  
##  [3] "No blast hit"                                                                                                                       
##  [4] "k__Fungi; p__Glomeromycota; c__Glomeromycetes; o__Glomerales; f__Glomeraceae; g__Glomus; s__sp"                                     
##  [5] "k__Fungi; p__Glomeromycota; c__Glomeromycetes; o__Diversisporales; f__Acaulosporaceae; g__Acaulospora; s__Acau16"                   
##  [6] "k__Fungi; p__Glomeromycota; c__Glomeromycetes; o__Glomerales; f__Glomeraceae; g__Glomus; s__MO_G47"                                 
##  [7] "k__Fungi; p__Glomeromycota; c__Glomeromycetes; o__Paraglomerales; f__Paraglomeraceae; g__Paraglomus; s__laccatum"                   
##  [8] "k__Fungi; p__Glomeromycota; c__Glomeromycetes; o__Glomerales; f__Glomeraceae; g__Glomus; s__acnaGlo2"                               
##  [9] "k__Fungi; p__Glomeromycota; c__Glomeromycetes; o__Archaeosporales; f__Archaeosporaceae; g__Archaeospora; s__MO_Ar1"                 
## [10] "k__Fungi; p__Glomeromycota; c__Glomeromycetes; o__Archaeosporales; f__Archaeosporaceae; g__Archaeospora; s__sp"                     
## [11] "k__Fungi; p__Glomeromycota; c__Glomeromycetes; o__Archaeosporales; f__Ambisporaceae; g__Ambispora; s__leptoticha"                   
## [12] "k__Fungi; p__Glomeromycota; c__Glomeromycetes; o__Glomerales; f__Glomeraceae; g__Glomus; s__Whitfield_type_7"                       
## [13] "k__Fungi; p__Glomeromycota; c__Glomeromycetes; o__Glomerales; f__Claroideoglomeraceae; g__Claroideoglomus; s__Torrecillas12b_Glo_G5"
## [14] "k__Fungi; p__Glomeromycota; c__Glomeromycetes; o__Glomerales; f__Glomeraceae; g__Glomus; s__Sanchez_Castro12b_GLO12"                
## [15] "k__Fungi; p__Glomeromycota; c__Glomeromycetes; o__Glomerales; f__Glomeraceae; g__Glomus; s__Winther07_D"                            
## [16] "k__Fungi; p__Glomeromycota; c__Glomeromycetes; o__Glomerales; f__Glomeraceae; g__Glomus; s__MO_G18"                                 
## [17] "k__Fungi; p__Glomeromycota; c__Glomeromycetes; o__Glomerales; f__Glomeraceae; g__Glomus; s__intraradices"                           
## [18] "k__Fungi; p__Glomeromycota; c__Glomeromycetes; o__Paraglomerales; f__Paraglomeraceae; g__Paraglomus; s__Alguacil12a_Para_1"         
## [19] "k__Fungi; p__Glomeromycota; c__Glomeromycetes; o__Glomerales; f__Claroideoglomeraceae; g__Claroideoglomus; s__Douhan9"              
## [20] "k__Fungi; p__Glomeromycota; c__Glomeromycetes; o__Glomerales; f__Claroideoglomeraceae; g__Claroideoglomus; s__Alguacil12b_GLO_G3"   
## [21] "k__Fungi; p__Glomeromycota; c__Glomeromycetes; o__Glomerales; f__Glomeraceae; g__Glomus; s__Glo39"                                  
## [22] "k__Fungi; p__Glomeromycota; c__Glomeromycetes; o__Glomerales; f__Glomeraceae; g__Glomus; s__Alguacil10_Glo1"                        
## [23] "k__Fungi; p__Glomeromycota; c__Glomeromycetes; o__Glomerales; f__Glomeraceae; g__Glomus; s__Ligrone07_sp"                           
## [24] "k__Fungi; p__Glomeromycota; c__Glomeromycetes; o__Diversisporales; f__Acaulosporaceae; g__Kuklospora; s__PT6"                       
## [25] "k__Fungi; p__Glomeromycota; c__Glomeromycetes; o__Diversisporales; f__Diversisporaceae; g__Diversispora; s__sp"                     
## [26] "k__Fungi; p__Glomeromycota; c__Glomeromycetes; o__Glomerales; f__Glomeraceae; g__Glomus; s__MO_G41"                                 
## [27] "k__Fungi; p__Glomeromycota; c__Glomeromycetes; o__Glomerales; f__Claroideoglomeraceae; g__Claroideoglomus; s__MO_GB1"               
## [28] "k__Fungi; p__Glomeromycota; c__Glomeromycetes; o__Glomerales; f__Glomeraceae; g__Glomus; s__A1"                                     
## [29] "k__Fungi; p__Glomeromycota; c__Glomeromycetes; o__Glomerales; f__Glomeraceae; g__Glomus; s__Liu2012b_Phylo_5"                       
## [30] "k__Fungi; p__Glomeromycota; c__Glomeromycetes; o__Glomerales; f__Glomeraceae; g__Glomus; s__Yamato09_C1"                            
## [31] "k__Fungi; p__Glomeromycota; c__Glomeromycetes; o__Glomerales; f__Claroideoglomeraceae; g__Claroideoglomus; s__acnaGlo7"             
## [32] "k__Fungi; p__Glomeromycota; c__Glomeromycetes; o__Glomerales; f__Glomeraceae; g__Glomus; s__MO_G7"                                  
## [33] "k__Fungi; p__Glomeromycota; c__Glomeromycetes; o__Glomerales; f__Glomeraceae; g__Glomus; s__Alguacil10_Glo6"                        
## [34] "k__Fungi; p__Glomeromycota; c__Glomeromycetes; o__Paraglomerales; f__Paraglomeraceae; g__Paraglomus; s__Alguacil12b_PARA1"          
## [35] "k__Fungi; p__Glomeromycota; c__Glomeromycetes; o__Glomerales; f__Claroideoglomeraceae; g__Claroideoglomus; s__Glo58"                
## [36] "k__Fungi; p__Glomeromycota; c__Glomeromycetes; o__Glomerales; f__Claroideoglomeraceae; g__Claroideoglomus; s__lamellosum"           
## [37] "k__Fungi; p__Glomeromycota; c__Glomeromycetes; o__Glomerales; f__Glomeraceae; g__Glomus; s__VeGlo18"                                
## [38] "k__Fungi; p__Glomeromycota; c__Glomeromycetes; o__Archaeosporales; f__Ambisporaceae; g__Ambispora; s__Liu2012a_Ar_1"                
## [39] "k__Fungi; p__Glomeromycota; c__Glomeromycetes; o__Diversisporales; f__Diversisporaceae; g__Diversispora; s__Div"                    
## [40] "k__Fungi; p__Glomeromycota; c__Glomeromycetes; o__Archaeosporales; f__Ambisporaceae; g__Ambispora; s__fennica"                      
## [41] "k__Fungi; p__Glomeromycota; c__Glomeromycetes; o__Glomerales; f__Glomeraceae; g__Glomus; s__Voyriella_parviflora_symbiont_type_1"   
## [42] "k__Fungi; p__Glomeromycota; c__Glomeromycetes; o__Paraglomerales; f__Paraglomeraceae; g__Paraglomus; s__sp"                         
## [43] "k__Fungi; p__Glomeromycota; c__Glomeromycetes; o__Paraglomerales; f__Paraglomeraceae; g__Paraglomus; s__Para1_OTU2"                 
## [44] "k__Fungi; p__Glomeromycota; c__Glomeromycetes; o__Glomerales; f__Glomeraceae; g__Glomus; s__Whitfield_type_3"                       
## [45] "k__Fungi; p__Glomeromycota; c__Glomeromycetes; o__Archaeosporales; f__Archaeosporaceae; g__Archaeospora; s__Aca"                    
## [46] "k__Fungi; p__Glomeromycota; c__Glomeromycetes; o__Glomerales; f__Glomeraceae; g__Glomus; s__MO_G20"                                 
## [47] "k__Fungi; p__Glomeromycota; c__Glomeromycetes; o__Diversisporales; f__Acaulosporaceae; g__Acaulospora; s__sp"                       
## [48] "k__Fungi; p__Glomeromycota; c__Glomeromycetes; o__Glomerales; f__Glomeraceae; g__Glomus; s__Glo7"                                   
## [49] "k__Fungi; p__Glomeromycota; c__Glomeromycetes; o__Diversisporales; f__Diversisporaceae; g__Diversispora; s__MO_GC1"                 
## [50] "k__Fungi; p__Glomeromycota; c__Glomeromycetes; o__Glomerales; f__Glomeraceae; g__Glomus; s__Winther07_B"                            
## [51] "k__Fungi; p__Glomeromycota; c__Glomeromycetes; o__Glomerales; f__Glomeraceae; g__Glomus; s__mosseae"                                
## [52] "k__Fungi; p__Glomeromycota; c__Glomeromycetes; o__Glomerales; f__Glomeraceae; g__Glomus; s__caledonium"                             
## [53] "k__Fungi; p__Glomeromycota; c__Glomeromycetes; o__Paraglomerales; f__Paraglomeraceae; g__Paraglomus; s__brasilianum"                
## [54] "k__Fungi; p__Glomeromycota; c__Glomeromycetes; o__Glomerales; f__Glomeraceae; g__Glomus; s__Alguacil09b_Glo_G8"                     
## [55] "k__Fungi; p__Glomeromycota; c__Glomeromycetes; o__Archaeosporales; f__Archaeosporaceae; g__Archaeospora; s__Schechter08_Arch1"      
## [56] "k__Fungi; p__Glomeromycota; c__Glomeromycetes; o__Archaeosporales; f__Archaeosporaceae; g__Archaeospora; s__Wirsel_OTU21"           
## [57] "k__Fungi; p__Glomeromycota; c__Glomeromycetes; o__Glomerales; f__Claroideoglomeraceae; g__Claroideoglomus; s__Glo71"                
## [58] "k__Fungi; p__Glomeromycota; c__Glomeromycetes; o__Diversisporales; f__Diversisporaceae; g__Diversispora; s__spurca"                 
## [59] "k__Fungi; p__Glomeromycota; c__Glomeromycetes; o__Glomerales; f__Claroideoglomeraceae; g__Claroideoglomus; s__Glo59"                
## [60] "k__Fungi; p__Glomeromycota; c__Glomeromycetes; o__Glomerales; f__Glomeraceae; g__Glomus; s__Wirsel_OTU16"                           
## [61] "k__Fungi; p__Glomeromycota; c__Glomeromycetes; o__Glomerales; f__Glomeraceae; g__Glomus; s__Glo32"                                  
## [62] "k__Fungi; p__Glomeromycota; c__Glomeromycetes; o__Diversisporales; f__Acaulosporaceae; g__Acaulospora; s__MO_A8"                    
## [63] "k__Fungi; p__Glomeromycota; c__Glomeromycetes; o__Glomerales; f__Glomeraceae; g__Glomus; s__MO_G27"                                 
## [64] "k__Fungi; p__Glomeromycota; c__Glomeromycetes; o__Glomerales; f__Glomeraceae; g__Glomus; s__MO_G5"
drop <- c("No blast hit")
ssu <- ssu[-which(ssu$ssutaxonomy %in% drop), ]

## alternatively, can use filtering to drop these and avoid rewriting things
ssu <- ssu.save %>% filter(ssufamily != "Geosiphonaceae" & ssutaxonomy != "No blast hit")

ssu <- ssu[complete.cases(ssu), ]
unique(ssu$ssuspecies)
##  [1] "NES27"                               
##  [2] "MH3_B"                               
##  [3] "sp"                                  
##  [4] "Acau16"                              
##  [5] "MO_G47"                              
##  [6] "laccatum"                            
##  [7] "acnaGlo2"                            
##  [8] "MO_Ar1"                              
##  [9] "leptoticha"                          
## [10] "Whitfield_type_7"                    
## [11] "Torrecillas12b_Glo_G5"               
## [12] "Sanchez_Castro12b_GLO12"             
## [13] "Winther07_D"                         
## [14] "MO_G18"                              
## [15] "intraradices"                        
## [16] "Alguacil12a_Para_1"                  
## [17] "Douhan9"                             
## [18] "Alguacil12b_GLO_G3"                  
## [19] "Glo39"                               
## [20] "Alguacil10_Glo1"                     
## [21] "Ligrone07_sp"                        
## [22] "PT6"                                 
## [23] "MO_G41"                              
## [24] "MO_GB1"                              
## [25] "A1"                                  
## [26] "Liu2012b_Phylo_5"                    
## [27] "Yamato09_C1"                         
## [28] "acnaGlo7"                            
## [29] "MO_G7"                               
## [30] "Alguacil10_Glo6"                     
## [31] "Alguacil12b_PARA1"                   
## [32] "Glo58"                               
## [33] "lamellosum"                          
## [34] "VeGlo18"                             
## [35] "Liu2012a_Ar_1"                       
## [36] "Div"                                 
## [37] "fennica"                             
## [38] "Voyriella_parviflora_symbiont_type_1"
## [39] "Para1_OTU2"                          
## [40] "Whitfield_type_3"                    
## [41] "Aca"                                 
## [42] "MO_G20"                              
## [43] "Glo7"                                
## [44] "MO_GC1"                              
## [45] "Winther07_B"                         
## [46] "mosseae"                             
## [47] "caledonium"                          
## [48] "brasilianum"                         
## [49] "Alguacil09b_Glo_G8"                  
## [50] "Schechter08_Arch1"                   
## [51] "Wirsel_OTU21"                        
## [52] "Glo71"                               
## [53] "spurca"                              
## [54] "Glo59"                               
## [55] "Wirsel_OTU16"                        
## [56] "Glo32"                               
## [57] "MO_A8"                               
## [58] "MO_G27"                              
## [59] "MO_G5"

Add in functional groups

ssul <- melt(ssu, variable.name = "ID", value.name = "ssureads")
str(ssul)
## 'data.frame':    17568 obs. of  11 variables:
##  $ ssuotu     : chr  "denovo538" "denovo1813859" "denovo34" "denovo11952" ...
##  $ ssutaxonomy: chr  "k__Fungi; p__Glomeromycota; c__Glomeromycetes; o__Glomerales; f__Glomeraceae; g__Glomus; s__NES27" "k__Fungi; p__Glomeromycota; c__Glomeromycetes; o__Glomerales; f__Glomeraceae; g__Glomus; s__MH3_B" "k__Fungi; p__Glomeromycota; c__Glomeromycetes; o__Glomerales; f__Glomeraceae; g__Glomus; s__NES27" "k__Fungi; p__Glomeromycota; c__Glomeromycetes; o__Glomerales; f__Glomeraceae; g__Glomus; s__sp" ...
##  $ ssukingdom : chr  "Fungi" "Fungi" "Fungi" "Fungi" ...
##  $ ssuphylum  : chr  "Glomeromycota" "Glomeromycota" "Glomeromycota" "Glomeromycota" ...
##  $ ssuclass   : chr  "Glomeromycetes" "Glomeromycetes" "Glomeromycetes" "Glomeromycetes" ...
##  $ ssuorder   : chr  "Glomerales" "Glomerales" "Glomerales" "Glomerales" ...
##  $ ssufamily  : chr  "Glomeraceae" "Glomeraceae" "Glomeraceae" "Glomeraceae" ...
##  $ ssugenus   : chr  "Glomus" "Glomus" "Glomus" "Glomus" ...
##  $ ssuspecies : chr  "NES27" "MH3_B" "NES27" "sp" ...
##  $ ID         : Factor w/ 61 levels "MM933M","MM936M",..: 1 1 1 1 1 1 1 1 1 1 ...
##  $ ssureads   : num  5.29 4.33 0 0 3.78 ...

OTU richnesses and read abundance for each taxonomic group

Below is in format for graphing

ssugenusREADRICH <- data.frame(ssul %>% group_by(ID, ssugenus) %>% summarise(OTU_Richness_Sample = length(unique(ssuotu[ssureads > 
    0])), Read_Abundance_Sample = sum(ssureads)))

Make SSU data frame with genus level OTU richness and Read abundance for each genus

ssugenera <- data.frame(ssul %>% group_by(ID, ssugenus) %>% summarise(Genus_OTU_Richness = length(unique(ssuotu[ssureads > 
    0])), Genus_Read_Abundance = sum(ssureads)))

Add metadata into ssu genus

colnames(metassu)  # need to edit for appropriate headers
##  [1] "#SampleID"            "BarcodeSequence"      "LinkerPrimerSequence"
##  [4] "Location"             "Project"              "Year"                
##  [7] "Site"                 "Treat"                "Rep"                 
## [10] "HL"                   "TSDel"                "Description"         
## [13] "SampleID2"
names(metassu)[names(metassu) == "SampleID2"] <- "ID"

ssugenera <- ssugenera %>% left_join(metassu, by = "ID")
str(ssugenera)
## 'data.frame':    488 obs. of  16 variables:
##  $ ID                  : chr  "MM933M" "MM933M" "MM933M" "MM933M" ...
##  $ ssugenus            : chr  "Acaulospora" "Ambispora" "Archaeospora" "Claroideoglomus" ...
##  $ Genus_OTU_Richness  : int  2 5 17 15 2 65 2 22 1 5 ...
##  $ Genus_Read_Abundance: num  6.66 23.94 157.59 99.62 23.2 ...
##  $ #SampleID           : int  5 5 5 5 5 5 5 5 3 3 ...
##  $ BarcodeSequence     : chr  "ATCCCGTATCGATTGG" "ATCCCGTATCGATTGG" "ATCCCGTATCGATTGG" "ATCCCGTATCGATTGG" ...
##  $ LinkerPrimerSequence: logi  NA NA NA NA NA NA ...
##  $ Location            : chr  "E1" "E1" "E1" "E1" ...
##  $ Project             : chr  "MM-Salvage" "MM-Salvage" "MM-Salvage" "MM-Salvage" ...
##  $ Year                : chr  "2015" "2015" "2015" "2015" ...
##  $ Site                : chr  "WL" "WL" "WL" "WL" ...
##  $ Treat               : Factor w/ 6 levels "C","D","F","O",..: 4 4 4 4 4 4 4 4 4 4 ...
##  $ Rep                 : Factor w/ 6 levels "1","2","3","4",..: 5 5 5 5 5 5 5 5 3 3 ...
##  $ HL                  : num  0.264 0.264 0.264 0.264 0.264 ...
##  $ TSDel               : chr  NA NA NA NA ...
##  $ Description         : chr  "RecipientPre" "RecipientPre" "RecipientPre" "RecipientPre" ...

Adds meta data to data for box plot

ssugenusREADRICH <- ssugenusREADRICH %>% left_join(metassu, by = "ID")
str(ssugenusREADRICH)
## 'data.frame':    488 obs. of  16 variables:
##  $ ID                   : chr  "MM933M" "MM933M" "MM933M" "MM933M" ...
##  $ ssugenus             : chr  "Acaulospora" "Ambispora" "Archaeospora" "Claroideoglomus" ...
##  $ OTU_Richness_Sample  : int  2 5 17 15 2 65 2 22 1 5 ...
##  $ Read_Abundance_Sample: num  6.66 23.94 157.59 99.62 23.2 ...
##  $ #SampleID            : int  5 5 5 5 5 5 5 5 3 3 ...
##  $ BarcodeSequence      : chr  "ATCCCGTATCGATTGG" "ATCCCGTATCGATTGG" "ATCCCGTATCGATTGG" "ATCCCGTATCGATTGG" ...
##  $ LinkerPrimerSequence : logi  NA NA NA NA NA NA ...
##  $ Location             : chr  "E1" "E1" "E1" "E1" ...
##  $ Project              : chr  "MM-Salvage" "MM-Salvage" "MM-Salvage" "MM-Salvage" ...
##  $ Year                 : chr  "2015" "2015" "2015" "2015" ...
##  $ Site                 : chr  "WL" "WL" "WL" "WL" ...
##  $ Treat                : Factor w/ 6 levels "C","D","F","O",..: 4 4 4 4 4 4 4 4 4 4 ...
##  $ Rep                  : Factor w/ 6 levels "1","2","3","4",..: 5 5 5 5 5 5 5 5 3 3 ...
##  $ HL                   : num  0.264 0.264 0.264 0.264 0.264 ...
##  $ TSDel                : chr  NA NA NA NA ...
##  $ Description          : chr  "RecipientPre" "RecipientPre" "RecipientPre" "RecipientPre" ...
filename <- paste0(date, "_ssu_genera_read_rich.csv")
write.csv(ssugenusREADRICH, file = filename, row.names = FALSE)

move ID to single line and genus to the length

### For ssufamily, moves ID to single line and genus to the length
LssugenusOTURICH <- dcast(ssugenera, ID ~ ssugenus, value.var = "Genus_OTU_Richness", 
    fun.aggregate = sum)
LssugenusOTUREAD <- dcast(ssugenera, ID ~ ssugenus, value.var = "Genus_Read_Abundance", 
    fun.aggregate = sum)

head(LssugenusOTUREAD)
##         ID Acaulospora Ambispora Archaeospora Claroideoglomus Diversispora
## 1   MM933M      6.6618   23.9417     157.5890         99.6242      23.1980
## 2 MM9362CU     12.2907   33.2473      63.2683        118.7108      18.3924
## 3   MM936M      3.3956   25.9602      67.7783         91.5404      12.2380
## 4   MM941M     21.9872   14.6280      47.8637         93.6373      13.9450
## 5   MM943M     14.2272   23.4654      92.1497         74.7298      22.9824
## 6   MMS700     11.9400   20.1172      90.3423        106.1814      18.6176
##     Glomus Kuklospora Paraglomus
## 1 387.3753     5.7632   129.6337
## 2 464.2526     3.1134   172.7572
## 3 406.1616     6.7912    81.1362
## 4 503.8588    11.0910    81.3654
## 5 317.0982    30.6556   103.0184
## 6 429.0388    10.8650   118.3515

Genera reads and richness together ###WORKED ON THIS 20180826

REPLACED Functional.Group with ssugenus

Work on changing to species in the next iteration

ssugenusREADRICHlong <- data.frame(ssul %>% group_by(ID) %>% summarise(OTU_Richness_Sample = length(unique(ssuotu[ssureads > 
    0])), Acaulospora_Richness = length(unique(ssuotu[ssureads > 0 & ssugenus == 
    "Acaulospora"])), Ambispora_Richness = length(unique(ssuotu[ssureads > 0 & 
    ssugenus == "Ambispora"])), Archaeospora_Richness = length(unique(ssuotu[ssureads > 
    0 & ssugenus == "Archaeospora"])), Claroideoglomus_Richness = length(unique(ssuotu[ssureads > 
    0 & ssugenus == "Claroideoglomus"])), Diversispora_Richness = length(unique(ssuotu[ssureads > 
    0 & ssugenus == "Diversispora"])), Glomus_Richness = length(unique(ssuotu[ssureads > 
    0 & ssugenus == "Glomus"])), Kuklospora_Richness = length(unique(ssuotu[ssureads > 
    0 & ssugenus == "Kuklospora"])), Paraglomus_Richness = length(unique(ssuotu[ssureads > 
    0 & ssugenus == "Paraglomus"])), Acaulospora_Read_Abundance = sum(ssureads[ssugenus == 
    "Acaulospora"]), Ambispora_Read_Abundance = sum(ssureads[ssugenus == "Ambispora"]), 
    Archaeospora_Read_Abundance = sum(ssureads[ssugenus == "Archaeospora"]), 
    Claroideoglomus_Read_Abundance = sum(ssureads[ssugenus == "Claroideoglomus"]), 
    Diversispora_Read_Abundance = sum(ssureads[ssugenus == "Diversispora"]), 
    Glomus_Read_Abundance = sum(ssureads[ssugenus == "Glomus"]), Kuklospora_Read_Abundance = sum(ssureads[ssugenus == 
        "Kuklospora"]), Paraglomus_Read_Abundance = sum(ssureads[ssugenus == 
        "Paraglomus"]), Read_Abundance_Sample = sum(ssureads)))
## View(ssugenusREADRICHlong)
ssugenusREADRICHlong <- merge(ssugenusREADRICHlong, metassu, by = "ID")

filename <- paste0(date, "_ssu_genus_read_rich.csv")
write.csv(ssugenusREADRICHlong, file = filename, row.names = FALSE)

str(ssugenusREADRICHlong)
## 'data.frame':    61 obs. of  31 variables:
##  $ ID                            : Factor w/ 61 levels "MM933M","MM936M",..: 1 5 2 3 4 6 7 8 9 10 ...
##  $ OTU_Richness_Sample           : int  130 167 105 121 99 151 158 161 138 148 ...
##  $ Acaulospora_Richness          : int  2 3 1 4 3 2 4 2 2 3 ...
##  $ Ambispora_Richness            : int  5 7 5 3 3 5 5 7 8 6 ...
##  $ Archaeospora_Richness         : int  17 16 12 10 11 15 13 13 8 13 ...
##  $ Claroideoglomus_Richness      : int  15 20 13 14 11 19 18 21 20 19 ...
##  $ Diversispora_Richness         : int  2 3 1 1 3 3 2 4 2 1 ...
##  $ Glomus_Richness               : int  65 85 59 74 51 83 85 80 79 77 ...
##  $ Kuklospora_Richness           : int  2 2 2 2 3 1 0 2 1 1 ...
##  $ Paraglomus_Richness           : int  22 31 12 13 14 23 31 32 18 28 ...
##  $ Acaulospora_Read_Abundance    : num  6.66 12.29 3.4 21.99 14.23 ...
##  $ Ambispora_Read_Abundance      : num  23.9 33.2 26 14.6 23.5 ...
##  $ Archaeospora_Read_Abundance   : num  157.6 63.3 67.8 47.9 92.1 ...
##  $ Claroideoglomus_Read_Abundance: num  99.6 118.7 91.5 93.6 74.7 ...
##  $ Diversispora_Read_Abundance   : num  23.2 18.4 12.2 13.9 23 ...
##  $ Glomus_Read_Abundance         : num  387 464 406 504 317 ...
##  $ Kuklospora_Read_Abundance     : num  5.76 3.11 6.79 11.09 30.66 ...
##  $ Paraglomus_Read_Abundance     : num  129.6 172.8 81.1 81.4 103 ...
##  $ Read_Abundance_Sample         : num  834 886 695 788 678 ...
##  $ #SampleID                     : int  5 4 3 1 2 47 48 49 50 51 ...
##  $ BarcodeSequence               : chr  "ATCCCGTATCGATTGG" "GCAACCTTTCGATTGG" "GCAACCTTAGAGTGTG" "GCAACCTTTGGGTGAT" ...
##  $ LinkerPrimerSequence          : logi  NA NA NA NA NA NA ...
##  $ Location                      : chr  "E1" "D1" "C1" "A1" ...
##  $ Project                       : chr  "MM-Salvage" "MM-Salvage" "MM-Salvage" "MM-Salvage" ...
##  $ Year                          : chr  "2015" "2015" "2015" "2015" ...
##  $ Site                          : chr  "WL" "WL" "WL" "WL" ...
##  $ Treat                         : Factor w/ 6 levels "C","D","F","O",..: 4 4 4 4 4 2 2 2 2 2 ...
##  $ Rep                           : Factor w/ 6 levels "1","2","3","4",..: 5 4 3 1 2 1 2 3 1 2 ...
##  $ HL                            : num  0.264 0.739 0.264 0.264 0.739 ...
##  $ TSDel                         : chr  NA NA NA NA ...
##  $ Description                   : chr  "RecipientPre" "RecipientPre" "RecipientPre" "RecipientPre" ...
## View(ssugenusREADRICHlong)

Non Unifrac distances, try with bray

WORK ON THIS

Permanova

Permanova tests for differences in composition among groups
Reminder - permanova is always based on pairwise distances/dissimilarities.
Can either nest Treat in Year with Year/Treat Or, can restrict permutations within year, using strata=‘Year’ Or, can remove Year, and focus on Treat and other grouping varuables, like Site

Also, dist.j<-vegdist(ssugenusREADRICHlong[,-1], method=‘jaccard’, na.rm=TRUE)

USE BRAY #dist.bc<-vegdist(ssugenusREADRICHlong[,-1], method=‘bray’, na.rm=TRUE)

Part 3 - multivariate community analyses and data visualization

library(dplyr)
library(reshape2)
library(tidyr)
library(vegan)
library(ggplot2)  # ggplot resource -http://rpubs.com/collnell/ggplot2
library(tidyverse)  # useful packages in 1 - dplyr, ggplot2, tidyr +

community matrix at order level with mapping data

otu_map<-read.csv(‘/Users/maltz/Desktop/RdataBreathe/otus_by_class.csv’)

head(ssugenusREADRICHlong)
##         ID OTU_Richness_Sample Acaulospora_Richness Ambispora_Richness
## 1   MM933M                 130                    2                  5
## 2 MM9362CU                 167                    3                  7
## 3   MM936M                 105                    1                  5
## 4   MM941M                 121                    4                  3
## 5   MM943M                  99                    3                  3
## 6   MMS700                 151                    2                  5
##   Archaeospora_Richness Claroideoglomus_Richness Diversispora_Richness
## 1                    17                       15                     2
## 2                    16                       20                     3
## 3                    12                       13                     1
## 4                    10                       14                     1
## 5                    11                       11                     3
## 6                    15                       19                     3
##   Glomus_Richness Kuklospora_Richness Paraglomus_Richness
## 1              65                   2                  22
## 2              85                   2                  31
## 3              59                   2                  12
## 4              74                   2                  13
## 5              51                   3                  14
## 6              83                   1                  23
##   Acaulospora_Read_Abundance Ambispora_Read_Abundance
## 1                     6.6618                  23.9417
## 2                    12.2907                  33.2473
## 3                     3.3956                  25.9602
## 4                    21.9872                  14.6280
## 5                    14.2272                  23.4654
## 6                    11.9400                  20.1172
##   Archaeospora_Read_Abundance Claroideoglomus_Read_Abundance
## 1                    157.5890                        99.6242
## 2                     63.2683                       118.7108
## 3                     67.7783                        91.5404
## 4                     47.8637                        93.6373
## 5                     92.1497                        74.7298
## 6                     90.3423                       106.1814
##   Diversispora_Read_Abundance Glomus_Read_Abundance
## 1                     23.1980              387.3753
## 2                     18.3924              464.2526
## 3                     12.2380              406.1616
## 4                     13.9450              503.8588
## 5                     22.9824              317.0982
## 6                     18.6176              429.0388
##   Kuklospora_Read_Abundance Paraglomus_Read_Abundance
## 1                    5.7632                  129.6337
## 2                    3.1134                  172.7572
## 3                    6.7912                   81.1362
## 4                   11.0910                   81.3654
## 5                   30.6556                  103.0184
## 6                   10.8650                  118.3515
##   Read_Abundance_Sample #SampleID  BarcodeSequence LinkerPrimerSequence
## 1              833.7869         5 ATCCCGTATCGATTGG                   NA
## 2              886.0327         4 GCAACCTTTCGATTGG                   NA
## 3              695.0015         3 GCAACCTTAGAGTGTG                   NA
## 4              788.3764         1 GCAACCTTTGGGTGAT                   NA
## 5              678.3267         2 GCAACCTTTACTGTGC                   NA
## 6              805.4538        47 GAAGATCCCTCTCAAG                   NA
##   Location    Project Year Site Treat Rep        HL TSDel   Description
## 1       E1 MM-Salvage 2015   WL     O   5 0.2644545  <NA>  RecipientPre
## 2       D1 MM-Salvage 2015   WL     O   4 0.7392792  <NA>  RecipientPre
## 3       C1 MM-Salvage 2015   WL     O   3 0.2644545  <NA>  RecipientPre
## 4       A1 MM-Salvage 2015   WL     O   1 0.2644545  <NA>  RecipientPre
## 5       B1 MM-Salvage 2015   WL     O   2 0.7392792  <NA>  RecipientPre
## 6       G6 MM-Salvage 2017   WL     D   1 0.4678646 B1A2C RecipientPost
colnames(ssugenusREADRICHlong)
##  [1] "ID"                             "OTU_Richness_Sample"           
##  [3] "Acaulospora_Richness"           "Ambispora_Richness"            
##  [5] "Archaeospora_Richness"          "Claroideoglomus_Richness"      
##  [7] "Diversispora_Richness"          "Glomus_Richness"               
##  [9] "Kuklospora_Richness"            "Paraglomus_Richness"           
## [11] "Acaulospora_Read_Abundance"     "Ambispora_Read_Abundance"      
## [13] "Archaeospora_Read_Abundance"    "Claroideoglomus_Read_Abundance"
## [15] "Diversispora_Read_Abundance"    "Glomus_Read_Abundance"         
## [17] "Kuklospora_Read_Abundance"      "Paraglomus_Read_Abundance"     
## [19] "Read_Abundance_Sample"          "#SampleID"                     
## [21] "BarcodeSequence"                "LinkerPrimerSequence"          
## [23] "Location"                       "Project"                       
## [25] "Year"                           "Site"                          
## [27] "Treat"                          "Rep"                           
## [29] "HL"                             "TSDel"                         
## [31] "Description"

Create the comm.grps from ssugenusREADRICHlong

To make the comm.grps I want the first column (col 1), then I don’t want (want to de-select) cols 2 through 22, and then I want the rest of the data = cols 23-30

This didn’t work:

comm.grps<-ssugenusREADRICHlong%>%dplyr::select(‘#SampleID’:Description) #mapping data #heads[1,-c(12:17)] ##This will not work because it uses numeric syntax and the text, only use one or the other (only select; or only numeric -c) comm.grps<-ssugenusREADRICHlong%>%dplyr::select[‘ID’:Description,-c(2:22)] #mapping data

#So I used the longhand way of typing out all the column names that I wanted to keep

`?`(dplyr::select)

comm.grps <- ssugenusREADRICHlong %>% dplyr::select("ID", "Location", "Year", 
    "Site", "Treat", "Rep", "TSDel", "Description")  #mapping data
colnames(comm.grps)
## [1] "ID"          "Location"    "Year"        "Site"        "Treat"      
## [6] "Rep"         "TSDel"       "Description"

WORK ON THIS!!!

Make the comm.mat ##Make the changes for this for species for the next iteration

comm.mat <- ssugenusREADRICHlong %>% dplyr::select(Acaulospora_Richness:Paraglomus_Richness)  # community matrix - all but mapping data

comparing ecological communities

diversity vs composition

abundance and richness are univariate response variables used to quantify communities

in multivariate analyses we have these variables for multiple entities

similarly, multivariate analyses have counterparts in univariate stats - t-test, ANOVA, mutliple regression

univariate analyses of diversity

list.files()  #shows what is in folder
##  [1] "2018Sep19_ssu_funcguild_read_rich.csv"               
##  [2] "2018Sep19_ssu_genera_read_rich.csv"                  
##  [3] "2018Sep19_ssu_genus_read_rich.csv"                   
##  [4] "2018Sep19_ssu_species_read_rich.csv"                 
##  [5] "2018Sep24_ssu_funcguild_read_rich.csv"               
##  [6] "2018Sep24_ssu_genera_read_rich.csv"                  
##  [7] "2018Sep24_ssu_species_read_rich.csv"                 
##  [8] "2018Sep27_ssu_funcguild_read_rich.csv"               
##  [9] "2018Sep27_ssu_genera_read_rich.csv"                  
## [10] "2018Sep27_ssu_genus_read_rich.csv"                   
## [11] "2018Sep27_ssu_species_read_rich.csv"                 
## [12] "2018Sep28_ssu_genera_read_rich.csv"                  
## [13] "2018Sep28_ssu_genus_read_rich.csv"                   
## [14] "AMF.Rmd"                                             
## [15] "AMF2.Rmd"                                            
## [16] "AMF4.Rmd"                                            
## [17] "AMF5.1.Rmd"                                          
## [18] "AMFModeling.Rmd"                                     
## [19] "AMFModelingRevise180913.Rmd"                         
## [20] "AMFModelingSpecies_1.Rmd"                            
## [21] "AMFModelingTaxa_3.Rmd"                               
## [22] "AMFModelingTaxa1 2.Rmd"                              
## [23] "AMFModelingTaxa1.Rmd"                                
## [24] "AMFModelingTaxa2.Rmd"                                
## [25] "AMFSubsetRevise180913.Rmd"                           
## [26] "bray_curtis_nmds_converted.txt"                      
## [27] "CN.csv"                                              
## [28] "CSS_table_sorted.txt"                                
## [29] "HL.xlsx"                                             
## [30] "HyphalSlidesDataInput_LD.xlsx"                       
## [31] "IdentifyingIssues"                                   
## [32] "map_SSU_Salvage_qiimeformat.txt"                     
## [33] "mapSal.txt"                                          
## [34] "mapSal18924.txt"                                     
## [35] "mapSal18927.txt"                                     
## [36] "rsconnect"                                           
## [37] "salhyp.txt"                                          
## [38] "SalvageHyphalExtractionSlidesMM.calc_in process.xlsx"
## [39] "SalvageProjectDescription"                           
## [40] "SSU_species.rtf"                                     
## [41] "Test_cache"                                          
## [42] "Test_files"                                          
## [43] "Test.html"                                           
## [44] "Test.Rmd"                                            
## [45] "Test27_cache"                                        
## [46] "Test27_files"                                        
## [47] "Test27.html"                                         
## [48] "Test27.Rmd"                                          
## [49] "Test28_cache"                                        
## [50] "Test28.Rmd"                                          
## [51] "TestSpecies_cache"                                   
## [52] "TestSpecies_files"                                   
## [53] "TestSpecies.html"                                    
## [54] "TestSpecies.Rmd"                                     
## [55] "TestSpeciesCorr2.Rmd"                                
## [56] "TestSpeciesCorrected_cache"                          
## [57] "TestSpeciesCorrected_files"                          
## [58] "TestSpeciesCorrected.html"                           
## [59] "TestSpeciesCorrected.Rmd"                            
## [60] "TypesOfAMF"
# View(comm.mat)
str(comm.grps)
## 'data.frame':    61 obs. of  8 variables:
##  $ ID         : Factor w/ 61 levels "MM933M","MM936M",..: 1 5 2 3 4 6 7 8 9 10 ...
##  $ Location   : chr  "E1" "D1" "C1" "A1" ...
##  $ Year       : chr  "2015" "2015" "2015" "2015" ...
##  $ Site       : chr  "WL" "WL" "WL" "WL" ...
##  $ Treat      : Factor w/ 6 levels "C","D","F","O",..: 4 4 4 4 4 2 2 2 2 2 ...
##  $ Rep        : Factor w/ 6 levels "1","2","3","4",..: 5 4 3 1 2 1 2 3 1 2 ...
##  $ TSDel      : chr  NA NA NA NA ...
##  $ Description: chr  "RecipientPre" "RecipientPre" "RecipientPre" "RecipientPre" ...
head(comm.grps)
##         ID Location Year Site Treat Rep TSDel   Description
## 1   MM933M       E1 2015   WL     O   5  <NA>  RecipientPre
## 2 MM9362CU       D1 2015   WL     O   4  <NA>  RecipientPre
## 3   MM936M       C1 2015   WL     O   3  <NA>  RecipientPre
## 4   MM941M       A1 2015   WL     O   1  <NA>  RecipientPre
## 5   MM943M       B1 2015   WL     O   2  <NA>  RecipientPre
## 6   MMS700       G6 2017   WL     D   1 B1A2C RecipientPost
str(comm.mat)
## 'data.frame':    61 obs. of  8 variables:
##  $ Acaulospora_Richness    : int  2 3 1 4 3 2 4 2 2 3 ...
##  $ Ambispora_Richness      : int  5 7 5 3 3 5 5 7 8 6 ...
##  $ Archaeospora_Richness   : int  17 16 12 10 11 15 13 13 8 13 ...
##  $ Claroideoglomus_Richness: int  15 20 13 14 11 19 18 21 20 19 ...
##  $ Diversispora_Richness   : int  2 3 1 1 3 3 2 4 2 1 ...
##  $ Glomus_Richness         : int  65 85 59 74 51 83 85 80 79 77 ...
##  $ Kuklospora_Richness     : int  2 2 2 2 3 1 0 2 1 1 ...
##  $ Paraglomus_Richness     : int  22 31 12 13 14 23 31 32 18 28 ...
head(comm.mat)
##   Acaulospora_Richness Ambispora_Richness Archaeospora_Richness
## 1                    2                  5                    17
## 2                    3                  7                    16
## 3                    1                  5                    12
## 4                    4                  3                    10
## 5                    3                  3                    11
## 6                    2                  5                    15
##   Claroideoglomus_Richness Diversispora_Richness Glomus_Richness
## 1                       15                     2              65
## 2                       20                     3              85
## 3                       13                     1              59
## 4                       14                     1              74
## 5                       11                     3              51
## 6                       19                     3              83
##   Kuklospora_Richness Paraglomus_Richness
## 1                   2                  22
## 2                   2                  31
## 3                   2                  12
## 4                   2                  13
## 5                   3                  14
## 6                   1                  23

does diversity vary across groups?

compute diversity indices

indices <- comm.grps

indices\(richness <- rowSums(comm.mat >0) #indices\)shannon <- diversity(comm.mat, index=‘shannon’)

indices$rarified <- c(rarefy(comm.mat, sample=50)) # rarefied diversity for a given sample size

indices <- comm.grps
indices$richness <- rowSums(comm.mat > 0)
indices$shannon <- diversity(comm.mat, index = "shannon")
View(indices)
# indices$rarified <- c(rarefy(comm.mat, sample=18103)) # rarefied diversity
# for a given sample size

WORK ON THIS!!!

Depending on what I put first in the model (Treat or Description) The terms are either significant (p-value) or not-significant.

Subset the data

prevpost.mat = comm.mat[comm.grps$Site == "WL", ]
prevpost.grps = comm.grps[comm.grps$Site == "WL", ]
distWL.bc <- vegdist(prevpost.mat, method = "bray", na.rm = TRUE)
WLTreat.div <- adonis2(distWL.bc ~ Description, strata = Treat, data = prevpost.grps, 
    permutations = 999, method = "bray")
WLTreat.div
## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
## 
## adonis2(formula = distWL.bc ~ Description, data = prevpost.grps, permutations = 999, method = "bray", strata = Treat)
##             Df SumOfSqs      R2      F Pr(>F)   
## Description  1 0.073503 0.35663 9.9779  0.002 **
## Residual    18 0.132599 0.64337                 
## Total       19 0.206101 1.00000                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
dvpostWL.mat = comm.mat[(comm.grps$Site == "WL" & comm.grps$Description == "RecipientPost") | 
    comm.grps$Site == "DS", ]
dvpostWL.grps = comm.grps[(comm.grps$Site == "WL" & comm.grps$Description == 
    "RecipientPost") | comm.grps$Site == "DS", ]
distDWL.bc <- vegdist(dvpostWL.mat, method = "bray", na.rm = TRUE)
DWLTreat.div <- adonis2(distDWL.bc ~ Description, strata = Treat, data = dvpostWL.grps, 
    permutations = 999, method = "bray")
DWLTreat.div
## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
## 
## adonis2(formula = distDWL.bc ~ Description, data = dvpostWL.grps, permutations = 999, method = "bray", strata = Treat)
##             Df SumOfSqs      R2      F Pr(>F)  
## Description  1 0.015150 0.19184 4.0355  0.025 *
## Residual    17 0.063823 0.80816                
## Total       18 0.078974 1.00000                
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
prevpost.mat = comm.mat[comm.grps$Site == "WL", ]
prevpost.grps = comm.grps[comm.grps$Site == "WL", ]
distWL.bc <- vegdist(prevpost.mat, method = "bray", na.rm = TRUE)
WLTreat.div <- adonis2(distWL.bc ~ Description, strata = Treat, data = prevpost.grps, 
    permutations = 999, method = "bray")
WLTreat.div
## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
## 
## adonis2(formula = distWL.bc ~ Description, data = prevpost.grps, permutations = 999, method = "bray", strata = Treat)
##             Df SumOfSqs      R2      F Pr(>F)   
## Description  1 0.073503 0.35663 9.9779  0.004 **
## Residual    18 0.132599 0.64337                 
## Total       19 0.206101 1.00000                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
dvpostWL.mat = comm.mat[(comm.grps$Site == "WL" & comm.grps$Description == "RecipientPost") | 
    comm.grps$Site == "DS", ]
dvpostWL.grps = comm.grps[(comm.grps$Site == "WL" & comm.grps$Description == 
    "RecipientPost") | comm.grps$Site == "DS", ]
distDWL.bc <- vegdist(dvpostWL.mat, method = "bray", na.rm = TRUE)
DWLTreat.div <- adonis2(distDWL.bc ~ Description, strata = Treat, data = dvpostWL.grps, 
    permutations = 999, method = "bray")
DWLTreat.div
## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
## 
## adonis2(formula = distDWL.bc ~ Description, data = dvpostWL.grps, permutations = 999, method = "bray", strata = Treat)
##             Df SumOfSqs      R2      F Pr(>F)  
## Description  1 0.015150 0.19184 4.0355  0.016 *
## Residual    17 0.063823 0.80816                
## Total       18 0.078974 1.00000                
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
prevpostPS.mat = comm.mat[comm.grps$Site == "PS", ]
prevpostPS.grps = comm.grps[comm.grps$Site == "PS", ]
distPS.bc <- vegdist(prevpostPS.mat, method = "bray", na.rm = TRUE)
PSTreat.div <- adonis2(distPS.bc ~ Description, strata = Treat, data = prevpostPS.grps, 
    permutations = 999, method = "bray")
PSTreat.div
## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
## 
## adonis2(formula = distPS.bc ~ Description, data = prevpostPS.grps, permutations = 999, method = "bray", strata = Treat)
##             Df SumOfSqs      R2      F Pr(>F)  
## Description  1 0.012597 0.17068 3.4988  0.033 *
## Residual    17 0.061209 0.82932                
## Total       18 0.073806 1.00000                
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
dim(prevpostPS.grps)
## [1] 19  8
# View(prevpostPS.grps)


dvpostPS.mat = comm.mat[(comm.grps$Site == "PS" & comm.grps$Description == "RecipientPost") | 
    comm.grps$Site == "DS", ]
dvpostPS.grps = comm.grps[(comm.grps$Site == "PS" & comm.grps$Description == 
    "RecipientPost") | comm.grps$Site == "DS", ]
distDPS.bc <- vegdist(dvpostPS.mat, method = "bray", na.rm = TRUE)
DPSTreat.div <- adonis2(distDPS.bc ~ Description, strata = Treat, data = dvpostPS.grps, 
    permutations = 999, method = "bray")
DPSTreat.div
## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
## 
## adonis2(formula = distDPS.bc ~ Description, data = dvpostPS.grps, permutations = 999, method = "bray", strata = Treat)
##             Df SumOfSqs     R2     F Pr(>F)  
## Description  1 0.008275 0.1119 2.268  0.092 .
## Residual    18 0.065677 0.8881               
## Total       19 0.073953 1.0000               
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
dim(dvpostPS.grps)
## [1] 20  8
prevpostHH.mat = comm.mat[comm.grps$Site == "HH", ]
prevpostHH.grps = comm.grps[comm.grps$Site == "HH", ]
distHH.bc <- vegdist(prevpostHH.mat, method = "bray", na.rm = TRUE)
HHTreat.div <- adonis2(distHH.bc ~ Description, strata = Treat, data = prevpostHH.grps, 
    permutations = 999, method = "bray")
HHTreat.div
## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
## 
## adonis2(formula = distHH.bc ~ Description, data = prevpostHH.grps, permutations = 999, method = "bray", strata = Treat)
##             Df SumOfSqs      R2      F Pr(>F)
## Description  1 0.001707 0.02229 0.3419   0.69
## Residual    15 0.074889 0.97771              
## Total       16 0.076596 1.00000
dim(prevpostHH.grps)
## [1] 17  8
dvpostHH.mat = comm.mat[(comm.grps$Site == "HH" & comm.grps$Description == "RecipientPost") | 
    comm.grps$Site == "HH", ]
dvpostHH.grps = comm.grps[(comm.grps$Site == "HH" & comm.grps$Description == 
    "RecipientPost") | comm.grps$Site == "HH", ]
distDHH.bc <- vegdist(dvpostHH.mat, method = "bray", na.rm = TRUE)
DHHTreat.div <- adonis2(distDHH.bc ~ Description, strata = Treat, data = dvpostHH.grps, 
    permutations = 999, method = "bray")
DHHTreat.div
## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
## 
## adonis2(formula = distDHH.bc ~ Description, data = dvpostHH.grps, permutations = 999, method = "bray", strata = Treat)
##             Df SumOfSqs      R2      F Pr(>F)
## Description  1 0.001707 0.02229 0.3419  0.702
## Residual    15 0.074889 0.97771              
## Total       16 0.076596 1.00000
dim(dvpostHH.grps)
## [1] 17  8
prevpostPS.mat = comm.mat[comm.grps$Site == "PS", ]
prevpostPS.grps = comm.grps[comm.grps$Site == "PS", ]
distPS.bc <- vegdist(prevpostPS.mat, method = "bray", na.rm = TRUE)
PSTreat.div <- adonis2(distPS.bc ~ Description, strata = Treat, data = prevpostPS.grps, 
    permutations = 999, method = "bray")
PSTreat.div
## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
## 
## adonis2(formula = distPS.bc ~ Description, data = prevpostPS.grps, permutations = 999, method = "bray", strata = Treat)
##             Df SumOfSqs      R2      F Pr(>F)  
## Description  1 0.012597 0.17068 3.4988  0.025 *
## Residual    17 0.061209 0.82932                
## Total       18 0.073806 1.00000                
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
dim(prevpostPS.grps)
## [1] 19  8
View(prevpostPS.grps)


dvpostPS.mat = comm.mat[(comm.grps$Site == "PS" & comm.grps$Description == "RecipientPost") | 
    comm.grps$Site == "DS", ]
dvpostPS.grps = comm.grps[(comm.grps$Site == "PS" & comm.grps$Description == 
    "RecipientPost") | comm.grps$Site == "DS", ]
distDPS.bc <- vegdist(dvpostPS.mat, method = "bray", na.rm = TRUE)
DPSTreat.div <- adonis2(distDPS.bc ~ Description, strata = Treat, data = dvpostPS.grps, 
    permutations = 999, method = "bray")
DPSTreat.div
## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
## 
## adonis2(formula = distDPS.bc ~ Description, data = dvpostPS.grps, permutations = 999, method = "bray", strata = Treat)
##             Df SumOfSqs     R2     F Pr(>F)  
## Description  1 0.008275 0.1119 2.268  0.098 .
## Residual    18 0.065677 0.8881               
## Total       19 0.073953 1.0000               
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
dim(dvpostPS.grps)
## [1] 20  8
postTreatPS.mat = comm.mat[comm.grps$Site == "PS", ]
postTreatPS.grps = comm.grps[comm.grps$Site == "PS", ]
distTPS.bc <- vegdist(postTreatPS.mat, method = "bray", na.rm = TRUE)

View(postTreatPS.grps)

contr.mineCD <- function(...) cbind(c(-1, 1, 0, 0, 0, 0))
PSTreatCD <- adonis2(distTPS.bc ~ Treat, data = postTreatPS.grps, permutations = 999, 
    method = "bray", contr.unordered = "contr.mineCD")
PSTreatCD
## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
## 
## adonis2(formula = distTPS.bc ~ Treat, data = postTreatPS.grps, permutations = 999, method = "bray", contr.unordered = "contr.mineCD")
##          Df SumOfSqs      R2      F Pr(>F)  
## Treat     5 0.034883 0.47264 2.3302  0.023 *
## Residual 13 0.038923 0.52736                
## Total    18 0.073806 1.00000                
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
contr.mineCT <- function(...) cbind(c(-1, 0, 1, 0, 0, 0))
PSTreatCT <- adonis2(distTPS.bc ~ Treat, data = postTreatPS.grps, permutations = 999, 
    method = "bray", contr.unordered = "contr.mineCT")
PSTreatCT
## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
## 
## adonis2(formula = distTPS.bc ~ Treat, data = postTreatPS.grps, permutations = 999, method = "bray", contr.unordered = "contr.mineCT")
##          Df SumOfSqs      R2      F Pr(>F)  
## Treat     5 0.034883 0.47264 2.3302  0.028 *
## Residual 13 0.038923 0.52736                
## Total    18 0.073806 1.00000                
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
contr.mineOC <- function(...) cbind(c(-1, 0, 0, 0, 0, 1))
PSTreatOC <- adonis2(distTPS.bc ~ Treat, data = postTreatPS.grps, permutations = 999, 
    method = "bray", contr.unordered = "contr.mineOC")
PSTreatOC
## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
## 
## adonis2(formula = distTPS.bc ~ Treat, data = postTreatPS.grps, permutations = 999, method = "bray", contr.unordered = "contr.mineOC")
##          Df SumOfSqs      R2      F Pr(>F)  
## Treat     5 0.034883 0.47264 2.3302  0.023 *
## Residual 13 0.038923 0.52736                
## Total    18 0.073806 1.00000                
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
postTreatPS.div <- adonis2(distPS.bc ~ Description, strata = Treat, data = postTreatPS.grps, 
    permutations = 999, method = "bray")
postTreatPS.div
## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
## 
## adonis2(formula = distPS.bc ~ Description, data = postTreatPS.grps, permutations = 999, method = "bray", strata = Treat)
##             Df SumOfSqs      R2      F Pr(>F)  
## Description  1 0.012597 0.17068 3.4988  0.028 *
## Residual    17 0.061209 0.82932                
## Total       18 0.073806 1.00000                
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
dim(postTreatPS.grps)
## [1] 19  8
View(postTreatPS.grps)


dvpostPS.mat = comm.mat[(comm.grps$Site == "PS" & comm.grps$Description == "RecipientPost") | 
    comm.grps$Site == "DS", ]
dvpostPS.grps = comm.grps[(comm.grps$Site == "PS" & comm.grps$Description == 
    "RecipientPost") | comm.grps$Site == "DS", ]
distDPS.bc <- vegdist(dvpostPS.mat, method = "bray", na.rm = TRUE)
DPSTreat.div <- adonis2(distDPS.bc ~ Description, strata = Treat, data = dvpostPS.grps, 
    permutations = 999, method = "bray")
DPSTreat.div
## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
## 
## adonis2(formula = distDPS.bc ~ Description, data = dvpostPS.grps, permutations = 999, method = "bray", strata = Treat)
##             Df SumOfSqs     R2     F Pr(>F)  
## Description  1 0.008275 0.1119 2.268  0.087 .
## Residual    18 0.065677 0.8881               
## Total       19 0.073953 1.0000               
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
dim(dvpostPS.grps)
## [1] 20  8
set.seed(304)
`?`(vegdist)

## are the results the same with other (non evolutionary) dissimiarlity
## indices? dist.j<-vegdist(ssugenusREADRICHlong[,-1], method='jaccard')
## dist.bc<-vegdist(ssugenusREADRICHlong[,-1], method='bray', na.rm=TRUE)
## dist.j<-vegdist(comm.mat, method='jaccard')
dist.bc <- vegdist(comm.mat, method = "bray", na.rm = TRUE)
AMTreat.div <- adonis2(dist.bc ~ Treat, data = comm.grps, permutations = 999, 
    method = "bray")
AMTreat.div
## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
## 
## adonis2(formula = dist.bc ~ Treat, data = comm.grps, permutations = 999, method = "bray")
##          Df SumOfSqs      R2     F Pr(>F)    
## Treat     5  0.09754 0.24065 3.486  0.001 ***
## Residual 55  0.30779 0.75935                 
## Total    60  0.40533 1.00000                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
AMTD.div <- adonis2(dist.bc ~ Treat + Description, data = comm.grps, permutations = 999, 
    method = "bray")
AMTD.div
## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
## 
## adonis2(formula = dist.bc ~ Treat + Description, data = comm.grps, permutations = 999, method = "bray")
##             Df SumOfSqs      R2      F Pr(>F)   
## Treat        5  0.09754 0.24065 3.5118  0.003 **
## Description  1  0.00782 0.01928 1.4071  0.227   
## Residual    54  0.29997 0.74007                 
## Total       60  0.40533 1.00000                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
AMTSD.div <- adonis2(dist.bc ~ Treat + Site * Description, data = comm.grps, 
    permutations = 999, method = "bray")
AMTSD.div
## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
## 
## adonis2(formula = dist.bc ~ Treat + Site * Description, data = comm.grps, permutations = 999, method = "bray")
##                  Df SumOfSqs      R2      F Pr(>F)    
## Treat             5  0.09754 0.24065 3.9234  0.001 ***
## Site              3  0.03184 0.07857 2.1348  0.057 .  
## Site:Description  2  0.02733 0.06742 2.7480  0.052 .  
## Residual         50  0.24861 0.61337                  
## Total            60  0.40533 1.00000                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# AMD_ST_J.div<-adonis2(dist.bc~Description+Site*Treat, data=comm.grps,
# permutations = 999, method='bray') AMD_ST_J.div

AMD_ST_BC.div <- adonis2(dist.bc ~ Description + Site + Treat, data = comm.grps, 
    permutations = 999, method = "bray")
AMD_ST_BC.div
## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
## 
## adonis2(formula = dist.bc ~ Description + Site + Treat, data = comm.grps, permutations = 999, method = "bray")
##             Df SumOfSqs      R2      F Pr(>F)   
## Description  2  0.07616 0.18790 7.1762  0.002 **
## Site         2  0.02444 0.06030 2.3030  0.065 . 
## Treat        4  0.02878 0.07101 1.3559  0.230   
## Residual    52  0.27594 0.68079                 
## Total       60  0.40533 1.00000                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Depending on what comes first in the model, we see different significance of these factors

AMTreat.div <- adonis2(dist.bc ~ Treat * Description, data = comm.grps, permutations = 999, 
    method = "bray")
AMTreat.div
## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
## 
## adonis2(formula = dist.bc ~ Treat * Description, data = comm.grps, permutations = 999, method = "bray")
##             Df SumOfSqs      R2      F Pr(>F)   
## Treat        5  0.09754 0.24065 3.5118  0.002 **
## Description  1  0.00782 0.01928 1.4071  0.231   
## Residual    54  0.29997 0.74007                 
## Total       60  0.40533 1.00000                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
amfMDS<-metaMDS(comm.mat, distance="bray", k=2, trymax=35, autotransform=TRUE) ##k is the number of dimensions
## Square root transformation
## Wisconsin double standardization
## Run 0 stress 0.1032977 
## Run 1 stress 0.1190699 
## Run 2 stress 0.1180244 
## Run 3 stress 0.129419 
## Run 4 stress 0.1026719 
## ... New best solution
## ... Procrustes: rmse 0.02425996  max resid 0.1500037 
## Run 5 stress 0.1144604 
## Run 6 stress 0.1221095 
## Run 7 stress 0.1212653 
## Run 8 stress 0.1283063 
## Run 9 stress 0.1080377 
## Run 10 stress 0.1182926 
## Run 11 stress 0.1220324 
## Run 12 stress 0.1245569 
## Run 13 stress 0.1218783 
## Run 14 stress 0.1067808 
## Run 15 stress 0.1033579 
## Run 16 stress 0.1062276 
## Run 17 stress 0.1033075 
## Run 18 stress 0.1358752 
## Run 19 stress 0.1033716 
## Run 20 stress 0.1257308 
## Run 21 stress 0.1144457 
## Run 22 stress 0.1040688 
## Run 23 stress 0.1342131 
## Run 24 stress 0.1206015 
## Run 25 stress 0.1087087 
## Run 26 stress 0.1164982 
## Run 27 stress 0.4024507 
## Run 28 stress 0.1224657 
## Run 29 stress 0.1147626 
## Run 30 stress 0.1223208 
## Run 31 stress 0.1304137 
## Run 32 stress 0.1166109 
## Run 33 stress 0.1023312 
## ... New best solution
## ... Procrustes: rmse 0.01416885  max resid 0.06962781 
## Run 34 stress 0.1036938 
## Run 35 stress 0.120806 
## *** No convergence -- monoMDS stopping criteria:
##     35: stress ratio > sratmax
amfMDS ##metaMDS takes eaither a distance matrix or your community matrix (then requires method for 'distance=')
## 
## Call:
## metaMDS(comm = comm.mat, distance = "bray", k = 2, trymax = 35,      autotransform = TRUE) 
## 
## global Multidimensional Scaling using monoMDS
## 
## Data:     wisconsin(sqrt(comm.mat)) 
## Distance: bray 
## 
## Dimensions: 2 
## Stress:     0.1023312 
## Stress type 1, weak ties
## No convergent solutions - best solution after 35 tries
## Scaling: centring, PC rotation, halfchange scaling 
## Species: expanded scores based on 'wisconsin(sqrt(comm.mat))'
stressplot(amfMDS)

#install.packages('ggplot2') ##plotting package
library(ggplot2)

NMDS1 <- amfMDS$points[,1] ##also found using scores(amfMDS)
NMDS2 <- amfMDS$points[,2]
?cbind
amf.plot<-cbind(comm.grps, NMDS1, NMDS2, comm.mat)
p<-ggplot(amf.plot, aes(NMDS1, NMDS2, color=Year))+
  geom_point(position=position_jitter(.1), shape=3)+##separates overlapping points
  stat_ellipse(type='t',size =1)+ ##draws 95% confidence interval ellipses
  theme_minimal()
p

plot<-ggplot(amf.plot, aes(NMDS1, NMDS2, color=Treat))+
  stat_ellipse(type='t',size =1)+
  theme_minimal()+geom_text(data=amf.plot,aes(NMDS1, NMDS2, label=Site), position=position_jitter(.35))+
  annotate("text", x=min(NMDS1), y=min(NMDS2), label=paste('Stress =',round(amfMDS$stress,3))) #add stress to plot
plot

Subset the community data matrix Number of samples that are left

summary(comm.grps)
##         ID       Location             Year               Site          
##  MM933M  : 1   Length:61          Length:61          Length:61         
##  MM936M  : 1   Class :character   Class :character   Class :character  
##  MM941M  : 1   Mode  :character   Mode  :character   Mode  :character  
##  MM943M  : 1                                                           
##  MM9362CU: 1                                                           
##  MMS700  : 1                                                           
##  (Other) :55                                                           
##  Treat  Rep       TSDel           Description       
##  C: 9   1:19   Length:61          Length:61         
##  D: 9   2:19   Class :character   Class :character  
##  F: 9   3:17   Mode  :character   Mode  :character  
##  O:17   4: 3                                        
##  S: 8   5: 2                                        
##  T: 9   6: 1                                        
## 
class(comm.mat)
## [1] "data.frame"
# subset the comm.mat and comm.grps x=comm.mat[comm.grps$Year == '2015',]
# y=comm.grps[comm.grps$Year == '2015',]

commented out some difficult lines

# what is x?

# SSMDS<-metaMDS(x, distance='bray', k = 2, trymax = 35, autotransform =
# TRUE) ##k is the number of dimensions SSMDS ##metaMDS takes eaither a
# distance matrix or your community matrix (then requires method for
# 'distance=')

# stressplot(SSMDS)

# install.packages('ggplot2') ##plotting package
library(ggplot2)

# NMDS1 <- SSMDS$points[,1] ##also found using scores(amfMDS) NMDS2 <-
# SSMDS$points[,2]
`?`(cbind)

The code was stopping at this point, and I created a new chunk below

SSamf.plot<-cbind(y, NMDS1, NMDS2, x)

Sp<-ggplot(SSamf.plot, aes(NMDS1, NMDS2, color=Site))+ geom_point(position=position_jitter(.1), shape=7)+##separates overlapping points # stat_ellipse(type=‘t’,size =1)+ ##draws 95% confidence interval ellipses theme_minimal() Sp

SSp<-ggplot(SSamf.plot, aes(NMDS1, NMDS2, color=Site))+ geom_point(position=position_jitter(.1), shape=3)+##separates overlapping points stat_ellipse(type=‘t’,size =1)+ ##draws 95% confidence interval ellipses theme_minimal() SSp

plot<-ggplot(SSamf.plot, aes(NMDS1, NMDS2, color=Treat))+ stat_ellipse(type=‘t’,size =1)+ theme_minimal()+geom_text(data=amf.plot,aes(NMDS1, NMDS2, label=Site), position=position_jitter(.35))+ annotate(“text”, x=min(NMDS1), y=min(NMDS2), label=paste(‘Stress =’,round(amfMDS$stress,3))) #add stress to plot plot

ad.bc<-adonis2(dist.bc~Treat+Metadata1+Metadata, data=grps, permutations=1000)

ad.bc

Examining the statistic support for increased rhizophilic richness in recipientpost plots

ssu.fg <- ssu %>% mutate(Functional.Group = case_when(ssufamily %in% c("Glomeraceae", 
    "Claroideoglomeraceae", "Paraglomeraceae") ~ "Rhizophilic", ssufamily %in% 
    c("Gigasporaceae", "Diversisporaceae") ~ "Edaphophilic", ssufamily %in% 
    c("Ambisporaceae", "Archaeosporaceae", "Acaulosporaceae") ~ "Ancestral"))

ssul <- melt(ssu.fg, variable.name = "ID", value.name = "ssureads")
str(ssul)
## 'data.frame':    17568 obs. of  12 variables:
##  $ ssuotu          : chr  "denovo538" "denovo1813859" "denovo34" "denovo11952" ...
##  $ ssutaxonomy     : chr  "k__Fungi; p__Glomeromycota; c__Glomeromycetes; o__Glomerales; f__Glomeraceae; g__Glomus; s__NES27" "k__Fungi; p__Glomeromycota; c__Glomeromycetes; o__Glomerales; f__Glomeraceae; g__Glomus; s__MH3_B" "k__Fungi; p__Glomeromycota; c__Glomeromycetes; o__Glomerales; f__Glomeraceae; g__Glomus; s__NES27" "k__Fungi; p__Glomeromycota; c__Glomeromycetes; o__Glomerales; f__Glomeraceae; g__Glomus; s__sp" ...
##  $ ssukingdom      : chr  "Fungi" "Fungi" "Fungi" "Fungi" ...
##  $ ssuphylum       : chr  "Glomeromycota" "Glomeromycota" "Glomeromycota" "Glomeromycota" ...
##  $ ssuclass        : chr  "Glomeromycetes" "Glomeromycetes" "Glomeromycetes" "Glomeromycetes" ...
##  $ ssuorder        : chr  "Glomerales" "Glomerales" "Glomerales" "Glomerales" ...
##  $ ssufamily       : chr  "Glomeraceae" "Glomeraceae" "Glomeraceae" "Glomeraceae" ...
##  $ ssugenus        : chr  "Glomus" "Glomus" "Glomus" "Glomus" ...
##  $ ssuspecies      : chr  "NES27" "MH3_B" "NES27" "sp" ...
##  $ Functional.Group: chr  "Rhizophilic" "Rhizophilic" "Rhizophilic" "Rhizophilic" ...
##  $ ID              : Factor w/ 61 levels "MM933M","MM936M",..: 1 1 1 1 1 1 1 1 1 1 ...
##  $ ssureads        : num  5.29 4.33 0 0 3.78 ...
ssuguildREADRICH <- data.frame(ssul %>% group_by(ID, Functional.Group) %>% summarise(OTU_Richness_Sample = length(unique(ssuotu[ssureads > 
    0])), Read_Abundance_Sample = sum(ssureads)))
ssufamily <- data.frame(ssul %>% group_by(ID, ssufamily) %>% summarise(Family_OTU_Richness = length(unique(ssuotu[ssureads > 
    0])), Family_Read_Abundance = sum(ssureads)))

colnames(metassu)  # need to edit for appropriate headers
##  [1] "#SampleID"            "BarcodeSequence"      "LinkerPrimerSequence"
##  [4] "Location"             "Project"              "Year"                
##  [7] "Site"                 "Treat"                "Rep"                 
## [10] "HL"                   "TSDel"                "Description"         
## [13] "ID"
names(metassu)[names(metassu) == "SampleID2"] <- "ID"

ssufamily <- ssufamily %>% left_join(metassu, by = "ID")
str(ssufamily)
## 'data.frame':    427 obs. of  16 variables:
##  $ ID                   : chr  "MM933M" "MM933M" "MM933M" "MM933M" ...
##  $ ssufamily            : chr  "Acaulosporaceae" "Ambisporaceae" "Archaeosporaceae" "Claroideoglomeraceae" ...
##  $ Family_OTU_Richness  : int  4 5 17 15 2 65 22 3 5 12 ...
##  $ Family_Read_Abundance: num  12.4 23.9 157.6 99.6 23.2 ...
##  $ #SampleID            : int  5 5 5 5 5 5 5 3 3 3 ...
##  $ BarcodeSequence      : chr  "ATCCCGTATCGATTGG" "ATCCCGTATCGATTGG" "ATCCCGTATCGATTGG" "ATCCCGTATCGATTGG" ...
##  $ LinkerPrimerSequence : logi  NA NA NA NA NA NA ...
##  $ Location             : chr  "E1" "E1" "E1" "E1" ...
##  $ Project              : chr  "MM-Salvage" "MM-Salvage" "MM-Salvage" "MM-Salvage" ...
##  $ Year                 : chr  "2015" "2015" "2015" "2015" ...
##  $ Site                 : chr  "WL" "WL" "WL" "WL" ...
##  $ Treat                : Factor w/ 6 levels "C","D","F","O",..: 4 4 4 4 4 4 4 4 4 4 ...
##  $ Rep                  : Factor w/ 6 levels "1","2","3","4",..: 5 5 5 5 5 5 5 3 3 3 ...
##  $ HL                   : num  0.264 0.264 0.264 0.264 0.264 ...
##  $ TSDel                : chr  NA NA NA NA ...
##  $ Description          : chr  "RecipientPre" "RecipientPre" "RecipientPre" "RecipientPre" ...
ssuguildREADRICH <- ssuguildREADRICH %>% left_join(metassu, by = "ID")
str(ssuguildREADRICH)
## 'data.frame':    183 obs. of  16 variables:
##  $ ID                   : chr  "MM933M" "MM933M" "MM933M" "MM936M" ...
##  $ Functional.Group     : chr  "Ancestral" "Edaphophilic" "Rhizophilic" "Ancestral" ...
##  $ OTU_Richness_Sample  : int  26 2 102 20 1 84 19 1 101 20 ...
##  $ Read_Abundance_Sample: num  194 23.2 616.6 103.9 12.2 ...
##  $ #SampleID            : int  5 5 5 3 3 3 1 1 1 2 ...
##  $ BarcodeSequence      : chr  "ATCCCGTATCGATTGG" "ATCCCGTATCGATTGG" "ATCCCGTATCGATTGG" "GCAACCTTAGAGTGTG" ...
##  $ LinkerPrimerSequence : logi  NA NA NA NA NA NA ...
##  $ Location             : chr  "E1" "E1" "E1" "C1" ...
##  $ Project              : chr  "MM-Salvage" "MM-Salvage" "MM-Salvage" "MM-Salvage" ...
##  $ Year                 : chr  "2015" "2015" "2015" "2015" ...
##  $ Site                 : chr  "WL" "WL" "WL" "WL" ...
##  $ Treat                : Factor w/ 6 levels "C","D","F","O",..: 4 4 4 4 4 4 4 4 4 4 ...
##  $ Rep                  : Factor w/ 6 levels "1","2","3","4",..: 5 5 5 3 3 3 1 1 1 2 ...
##  $ HL                   : num  0.264 0.264 0.264 0.264 0.264 ...
##  $ TSDel                : chr  NA NA NA NA ...
##  $ Description          : chr  "RecipientPre" "RecipientPre" "RecipientPre" "RecipientPre" ...
filename <- paste0(date, "_ssu_funcguild_read_rich.csv")
write.csv(ssuguildREADRICH, file = filename, row.names = FALSE)

LssufamilyOTURICH <- dcast(ssufamily, ID ~ ssufamily, value.var = "Family_OTU_Richness", 
    fun.aggregate = sum)
LssufamilyOTUREAD <- dcast(ssufamily, ID ~ ssufamily, value.var = "Family_Read_Abundance", 
    fun.aggregate = sum)

head(LssufamilyOTUREAD)
##         ID Acaulosporaceae Ambisporaceae Archaeosporaceae
## 1   MM933M         12.4250       23.9417         157.5890
## 2 MM9362CU         15.4041       33.2473          63.2683
## 3   MM936M         10.1868       25.9602          67.7783
## 4   MM941M         33.0782       14.6280          47.8637
## 5   MM943M         44.8828       23.4654          92.1497
## 6   MMS700         22.8050       20.1172          90.3423
##   Claroideoglomeraceae Diversisporaceae Glomeraceae Paraglomeraceae
## 1              99.6242          23.1980    387.3753        129.6337
## 2             118.7108          18.3924    464.2526        172.7572
## 3              91.5404          12.2380    406.1616         81.1362
## 4              93.6373          13.9450    503.8588         81.3654
## 5              74.7298          22.9824    317.0982        103.0184
## 6             106.1814          18.6176    429.0388        118.3515
ssuguildREADRICHlong <- data.frame(ssul %>% group_by(ID) %>% summarise(OTU_Richness_Sample = length(unique(ssuotu[ssureads > 
    0])), Ancestral_Richness = length(unique(ssuotu[ssureads > 0 & Functional.Group == 
    "Ancestral"])), Edaphophilic_Richness = length(unique(ssuotu[ssureads > 
    0 & Functional.Group == "Edaphophilic"])), Rhizophilic_Richness = length(unique(ssuotu[ssureads > 
    0 & Functional.Group == "Rhizophilic"])), Ancestral_Read_Abundance = sum(ssureads[Functional.Group == 
    "Ancestral"]), Edaphophilic_Read_Abundance = sum(ssureads[Functional.Group == 
    "Edaphophilic"]), Rhizophilic_Read_Abundance = sum(ssureads[Functional.Group == 
    "Rhizophilic"]), Read_Abundance_Sample = sum(ssureads)))
View(ssuguildREADRICHlong)
ssuguildREADRICHlong <- merge(ssuguildREADRICHlong, metassu, by = "ID")

filename <- paste0(date, "_ssu_funcguild_read_rich.csv")
write.csv(ssuguildREADRICHlong, file = filename, row.names = FALSE)

str(ssuguildREADRICHlong)
## 'data.frame':    61 obs. of  21 variables:
##  $ ID                         : Factor w/ 61 levels "MM933M","MM936M",..: 1 5 2 3 4 6 7 8 9 10 ...
##  $ OTU_Richness_Sample        : int  130 167 105 121 99 151 158 161 138 148 ...
##  $ Ancestral_Richness         : int  26 28 20 19 20 23 22 24 19 23 ...
##  $ Edaphophilic_Richness      : int  2 3 1 1 3 3 2 4 2 1 ...
##  $ Rhizophilic_Richness       : int  102 136 84 101 76 125 134 133 117 124 ...
##  $ Ancestral_Read_Abundance   : num  194 111.9 103.9 95.6 160.5 ...
##  $ Edaphophilic_Read_Abundance: num  23.2 18.4 12.2 13.9 23 ...
##  $ Rhizophilic_Read_Abundance : num  617 756 579 679 495 ...
##  $ Read_Abundance_Sample      : num  834 886 695 788 678 ...
##  $ #SampleID                  : int  5 4 3 1 2 47 48 49 50 51 ...
##  $ BarcodeSequence            : chr  "ATCCCGTATCGATTGG" "GCAACCTTTCGATTGG" "GCAACCTTAGAGTGTG" "GCAACCTTTGGGTGAT" ...
##  $ LinkerPrimerSequence       : logi  NA NA NA NA NA NA ...
##  $ Location                   : chr  "E1" "D1" "C1" "A1" ...
##  $ Project                    : chr  "MM-Salvage" "MM-Salvage" "MM-Salvage" "MM-Salvage" ...
##  $ Year                       : chr  "2015" "2015" "2015" "2015" ...
##  $ Site                       : chr  "WL" "WL" "WL" "WL" ...
##  $ Treat                      : Factor w/ 6 levels "C","D","F","O",..: 4 4 4 4 4 2 2 2 2 2 ...
##  $ Rep                        : Factor w/ 6 levels "1","2","3","4",..: 5 4 3 1 2 1 2 3 1 2 ...
##  $ HL                         : num  0.264 0.739 0.264 0.264 0.739 ...
##  $ TSDel                      : chr  NA NA NA NA ...
##  $ Description                : chr  "RecipientPre" "RecipientPre" "RecipientPre" "RecipientPre" ...
ssuguildREAD <- ssul %>% group_by(ID) %>% summarise(Ancestral_Read_Abundance = sum(ssureads[Functional.Group == 
    "Ancestral"]), Edaphophilic_Read_Abundance = sum(ssureads[Functional.Group == 
    "Edaphophilic"]), Rhizophilic_Read_Abundance = sum(ssureads[Functional.Group == 
    "Rhizophilic"]))


ssuguildRICH <- ssul %>% group_by(ID) %>% summarise(OTU_Richness_Sample = length(unique(ssuotu[ssureads > 
    0])), Ancestral_Richness = length(unique(ssuotu[ssureads > 0 & Functional.Group == 
    "Ancestral"])), Edaphophilic_Richness = length(unique(ssuotu[ssureads > 
    0 & Functional.Group == "Edaphophilic"])), Rhizophilic_Richness = length(unique(ssuotu[ssureads > 
    0 & Functional.Group == "Rhizophilic"])))
head(ssuguildRICH)
## # A tibble: 6 x 5
##   ID    OTU_Richness_Sa… Ancestral_Richn… Edaphophilic_Ri… Rhizophilic_Ri…
##   <fct>            <int>            <int>            <int>           <int>
## 1 MM93…              130               26                2             102
## 2 MM93…              105               20                1              84
## 3 MM94…              121               19                1             101
## 4 MM94…               99               20                3              76
## 5 MM93…              167               28                3             136
## 6 MMS7…              151               23                3             125
ssuguildRICH <- merge(ssuguildRICH, metassu, by = "ID")
ssuguildREAD <- merge(ssuguildREAD, metassu, by = "ID")
str(ssuguildRICH)
## 'data.frame':    61 obs. of  17 variables:
##  $ ID                   : Factor w/ 61 levels "MM933M","MM936M",..: 1 5 2 3 4 6 7 8 9 10 ...
##  $ OTU_Richness_Sample  : int  130 167 105 121 99 151 158 161 138 148 ...
##  $ Ancestral_Richness   : int  26 28 20 19 20 23 22 24 19 23 ...
##  $ Edaphophilic_Richness: int  2 3 1 1 3 3 2 4 2 1 ...
##  $ Rhizophilic_Richness : int  102 136 84 101 76 125 134 133 117 124 ...
##  $ #SampleID            : int  5 4 3 1 2 47 48 49 50 51 ...
##  $ BarcodeSequence      : chr  "ATCCCGTATCGATTGG" "GCAACCTTTCGATTGG" "GCAACCTTAGAGTGTG" "GCAACCTTTGGGTGAT" ...
##  $ LinkerPrimerSequence : logi  NA NA NA NA NA NA ...
##  $ Location             : chr  "E1" "D1" "C1" "A1" ...
##  $ Project              : chr  "MM-Salvage" "MM-Salvage" "MM-Salvage" "MM-Salvage" ...
##  $ Year                 : chr  "2015" "2015" "2015" "2015" ...
##  $ Site                 : chr  "WL" "WL" "WL" "WL" ...
##  $ Treat                : Factor w/ 6 levels "C","D","F","O",..: 4 4 4 4 4 2 2 2 2 2 ...
##  $ Rep                  : Factor w/ 6 levels "1","2","3","4",..: 5 4 3 1 2 1 2 3 1 2 ...
##  $ HL                   : num  0.264 0.739 0.264 0.264 0.739 ...
##  $ TSDel                : chr  NA NA NA NA ...
##  $ Description          : chr  "RecipientPre" "RecipientPre" "RecipientPre" "RecipientPre" ...
str(ssuguildREAD)
## 'data.frame':    61 obs. of  16 variables:
##  $ ID                         : Factor w/ 61 levels "MM933M","MM936M",..: 1 5 2 3 4 6 7 8 9 10 ...
##  $ Ancestral_Read_Abundance   : num  194 111.9 103.9 95.6 160.5 ...
##  $ Edaphophilic_Read_Abundance: num  23.2 18.4 12.2 13.9 23 ...
##  $ Rhizophilic_Read_Abundance : num  617 756 579 679 495 ...
##  $ #SampleID                  : int  5 4 3 1 2 47 48 49 50 51 ...
##  $ BarcodeSequence            : chr  "ATCCCGTATCGATTGG" "GCAACCTTTCGATTGG" "GCAACCTTAGAGTGTG" "GCAACCTTTGGGTGAT" ...
##  $ LinkerPrimerSequence       : logi  NA NA NA NA NA NA ...
##  $ Location                   : chr  "E1" "D1" "C1" "A1" ...
##  $ Project                    : chr  "MM-Salvage" "MM-Salvage" "MM-Salvage" "MM-Salvage" ...
##  $ Year                       : chr  "2015" "2015" "2015" "2015" ...
##  $ Site                       : chr  "WL" "WL" "WL" "WL" ...
##  $ Treat                      : Factor w/ 6 levels "C","D","F","O",..: 4 4 4 4 4 2 2 2 2 2 ...
##  $ Rep                        : Factor w/ 6 levels "1","2","3","4",..: 5 4 3 1 2 1 2 3 1 2 ...
##  $ HL                         : num  0.264 0.739 0.264 0.264 0.739 ...
##  $ TSDel                      : chr  NA NA NA NA ...
##  $ Description                : chr  "RecipientPre" "RecipientPre" "RecipientPre" "RecipientPre" ...
y <- ssuguildRICH$Rhizophilic_Richness
str(y)
##  int [1:61] 102 136 84 101 76 125 134 133 117 124 ...
D <- ssuguildRICH$Description
S <- ssuguildRICH$Site

length(S)
## [1] 61
length(D)
## [1] 61
`?`(lm)
# lm(y~D)

SRhizRich <- lm(y ~ S)
str(SRhizRich)
## List of 13
##  $ coefficients : Named num [1:4] 125.4 12.6 6.18 3.2
##   ..- attr(*, "names")= chr [1:4] "(Intercept)" "SHH" "SPS" "SWL"
##  $ residuals    : Named num [1:61] -26.6 7.4 -44.6 -27.6 -52.6 ...
##   ..- attr(*, "names")= chr [1:61] "1" "2" "3" "4" ...
##  $ effects      : Named num [1:61] -1030.1 29.7 11.9 -6.4 -44.7 ...
##   ..- attr(*, "names")= chr [1:61] "(Intercept)" "SHH" "SPS" "SWL" ...
##  $ rank         : int 4
##  $ fitted.values: Named num [1:61] 129 129 129 129 129 ...
##   ..- attr(*, "names")= chr [1:61] "1" "2" "3" "4" ...
##  $ assign       : int [1:4] 0 1 1 1
##  $ qr           :List of 5
##   ..$ qr   : num [1:61, 1:4] -7.81 0.128 0.128 0.128 0.128 ...
##   .. ..- attr(*, "dimnames")=List of 2
##   .. .. ..$ : chr [1:61] "1" "2" "3" "4" ...
##   .. .. ..$ : chr [1:4] "(Intercept)" "SHH" "SPS" "SWL"
##   .. ..- attr(*, "assign")= int [1:4] 0 1 1 1
##   .. ..- attr(*, "contrasts")=List of 1
##   .. .. ..$ S: chr "contr.treatment"
##   ..$ qraux: num [1:4] 1.13 1.07 1.11 1.07
##   ..$ pivot: int [1:4] 1 2 3 4
##   ..$ tol  : num 1e-07
##   ..$ rank : int 4
##   ..- attr(*, "class")= chr "qr"
##  $ df.residual  : int 57
##  $ contrasts    :List of 1
##   ..$ S: chr "contr.treatment"
##  $ xlevels      :List of 1
##   ..$ S: chr [1:4] "DS" "HH" "PS" "WL"
##  $ call         : language lm(formula = y ~ S)
##  $ terms        :Classes 'terms', 'formula'  language y ~ S
##   .. ..- attr(*, "variables")= language list(y, S)
##   .. ..- attr(*, "factors")= int [1:2, 1] 0 1
##   .. .. ..- attr(*, "dimnames")=List of 2
##   .. .. .. ..$ : chr [1:2] "y" "S"
##   .. .. .. ..$ : chr "S"
##   .. ..- attr(*, "term.labels")= chr "S"
##   .. ..- attr(*, "order")= int 1
##   .. ..- attr(*, "intercept")= int 1
##   .. ..- attr(*, "response")= int 1
##   .. ..- attr(*, ".Environment")=<environment: R_GlobalEnv> 
##   .. ..- attr(*, "predvars")= language list(y, S)
##   .. ..- attr(*, "dataClasses")= Named chr [1:2] "numeric" "character"
##   .. .. ..- attr(*, "names")= chr [1:2] "y" "S"
##  $ model        :'data.frame':   61 obs. of  2 variables:
##   ..$ y: int [1:61] 102 136 84 101 76 125 134 133 117 124 ...
##   ..$ S: chr [1:61] "WL" "WL" "WL" "WL" ...
##   ..- attr(*, "terms")=Classes 'terms', 'formula'  language y ~ S
##   .. .. ..- attr(*, "variables")= language list(y, S)
##   .. .. ..- attr(*, "factors")= int [1:2, 1] 0 1
##   .. .. .. ..- attr(*, "dimnames")=List of 2
##   .. .. .. .. ..$ : chr [1:2] "y" "S"
##   .. .. .. .. ..$ : chr "S"
##   .. .. ..- attr(*, "term.labels")= chr "S"
##   .. .. ..- attr(*, "order")= int 1
##   .. .. ..- attr(*, "intercept")= int 1
##   .. .. ..- attr(*, "response")= int 1
##   .. .. ..- attr(*, ".Environment")=<environment: R_GlobalEnv> 
##   .. .. ..- attr(*, "predvars")= language list(y, S)
##   .. .. ..- attr(*, "dataClasses")= Named chr [1:2] "numeric" "character"
##   .. .. .. ..- attr(*, "names")= chr [1:2] "y" "S"
##  - attr(*, "class")= chr "lm"
DRhizRich <- lm(y ~ D)
str(DRhizRich)
## List of 13
##  $ coefficients : Named num [1:3] 125.4 11.44 -8.98
##   ..- attr(*, "names")= chr [1:3] "(Intercept)" "DRecipientPost" "DRecipientPre"
##  $ residuals    : Named num [1:61] -14.4 19.6 -32.4 -15.4 -40.4 ...
##   ..- attr(*, "names")= chr [1:61] "1" "2" "3" "4" ...
##  $ effects      : Named num [1:61] -1030.1 62.3 16.9 -13.5 -38.5 ...
##   ..- attr(*, "names")= chr [1:61] "(Intercept)" "DRecipientPost" "DRecipientPre" "" ...
##  $ rank         : int 3
##  $ fitted.values: Named num [1:61] 116 116 116 116 116 ...
##   ..- attr(*, "names")= chr [1:61] "1" "2" "3" "4" ...
##  $ assign       : int [1:3] 0 1 1
##  $ qr           :List of 5
##   ..$ qr   : num [1:61, 1:3] -7.81 0.128 0.128 0.128 0.128 ...
##   .. ..- attr(*, "dimnames")=List of 2
##   .. .. ..$ : chr [1:61] "1" "2" "3" "4" ...
##   .. .. ..$ : chr [1:3] "(Intercept)" "DRecipientPost" "DRecipientPre"
##   .. ..- attr(*, "assign")= int [1:3] 0 1 1
##   .. ..- attr(*, "contrasts")=List of 1
##   .. .. ..$ D: chr "contr.treatment"
##   ..$ qraux: num [1:3] 1.13 1.18 1.12
##   ..$ pivot: int [1:3] 1 2 3
##   ..$ tol  : num 1e-07
##   ..$ rank : int 3
##   ..- attr(*, "class")= chr "qr"
##  $ df.residual  : int 58
##  $ contrasts    :List of 1
##   ..$ D: chr "contr.treatment"
##  $ xlevels      :List of 1
##   ..$ D: chr [1:3] "Donor" "RecipientPost" "RecipientPre"
##  $ call         : language lm(formula = y ~ D)
##  $ terms        :Classes 'terms', 'formula'  language y ~ D
##   .. ..- attr(*, "variables")= language list(y, D)
##   .. ..- attr(*, "factors")= int [1:2, 1] 0 1
##   .. .. ..- attr(*, "dimnames")=List of 2
##   .. .. .. ..$ : chr [1:2] "y" "D"
##   .. .. .. ..$ : chr "D"
##   .. ..- attr(*, "term.labels")= chr "D"
##   .. ..- attr(*, "order")= int 1
##   .. ..- attr(*, "intercept")= int 1
##   .. ..- attr(*, "response")= int 1
##   .. ..- attr(*, ".Environment")=<environment: R_GlobalEnv> 
##   .. ..- attr(*, "predvars")= language list(y, D)
##   .. ..- attr(*, "dataClasses")= Named chr [1:2] "numeric" "character"
##   .. .. ..- attr(*, "names")= chr [1:2] "y" "D"
##  $ model        :'data.frame':   61 obs. of  2 variables:
##   ..$ y: int [1:61] 102 136 84 101 76 125 134 133 117 124 ...
##   ..$ D: chr [1:61] "RecipientPre" "RecipientPre" "RecipientPre" "RecipientPre" ...
##   ..- attr(*, "terms")=Classes 'terms', 'formula'  language y ~ D
##   .. .. ..- attr(*, "variables")= language list(y, D)
##   .. .. ..- attr(*, "factors")= int [1:2, 1] 0 1
##   .. .. .. ..- attr(*, "dimnames")=List of 2
##   .. .. .. .. ..$ : chr [1:2] "y" "D"
##   .. .. .. .. ..$ : chr "D"
##   .. .. ..- attr(*, "term.labels")= chr "D"
##   .. .. ..- attr(*, "order")= int 1
##   .. .. ..- attr(*, "intercept")= int 1
##   .. .. ..- attr(*, "response")= int 1
##   .. .. ..- attr(*, ".Environment")=<environment: R_GlobalEnv> 
##   .. .. ..- attr(*, "predvars")= language list(y, D)
##   .. .. ..- attr(*, "dataClasses")= Named chr [1:2] "numeric" "character"
##   .. .. .. ..- attr(*, "names")= chr [1:2] "y" "D"
##  - attr(*, "class")= chr "lm"
DSRhizRich <- lm(y ~ S * D)
str(DSRhizRich)
## List of 13
##  $ coefficients : Named num [1:12] 125.4 8.6 -3.4 -18.6 31.1 ...
##   ..- attr(*, "names")= chr [1:12] "(Intercept)" "SHH" "SPS" "SWL" ...
##  $ residuals    : Named num [1:61] -4.83 29.17 -22.83 -5.83 -30.83 ...
##   ..- attr(*, "names")= chr [1:61] "1" "2" "3" "4" ...
##  $ effects      : Named num [1:61] -1030.1 29.7 11.9 -6.4 58.6 ...
##   ..- attr(*, "names")= chr [1:61] "(Intercept)" "SHH" "SPS" "SWL" ...
##  $ rank         : int 7
##  $ fitted.values: Named num [1:61] 107 107 107 107 107 ...
##   ..- attr(*, "names")= chr [1:61] "1" "2" "3" "4" ...
##  $ assign       : int [1:12] 0 1 1 1 2 2 3 3 3 3 ...
##  $ qr           :List of 5
##   ..$ qr   : num [1:61, 1:12] -7.81 0.128 0.128 0.128 0.128 ...
##   .. ..- attr(*, "dimnames")=List of 2
##   .. .. ..$ : chr [1:61] "1" "2" "3" "4" ...
##   .. .. ..$ : chr [1:12] "(Intercept)" "SHH" "SPS" "SWL" ...
##   .. ..- attr(*, "assign")= int [1:12] 0 1 1 1 2 2 3 3 3 3 ...
##   .. ..- attr(*, "contrasts")=List of 2
##   .. .. ..$ S: chr "contr.treatment"
##   .. .. ..$ D: chr "contr.treatment"
##   ..$ qraux: num [1:12] 1.13 1.07 1.11 1.07 1.16 ...
##   ..$ pivot: int [1:12] 1 2 3 4 5 7 8 6 9 10 ...
##   ..$ tol  : num 1e-07
##   ..$ rank : int 7
##   ..- attr(*, "class")= chr "qr"
##  $ df.residual  : int 54
##  $ contrasts    :List of 2
##   ..$ S: chr "contr.treatment"
##   ..$ D: chr "contr.treatment"
##  $ xlevels      :List of 2
##   ..$ S: chr [1:4] "DS" "HH" "PS" "WL"
##   ..$ D: chr [1:3] "Donor" "RecipientPost" "RecipientPre"
##  $ call         : language lm(formula = y ~ S * D)
##  $ terms        :Classes 'terms', 'formula'  language y ~ S * D
##   .. ..- attr(*, "variables")= language list(y, S, D)
##   .. ..- attr(*, "factors")= int [1:3, 1:3] 0 1 0 0 0 1 0 1 1
##   .. .. ..- attr(*, "dimnames")=List of 2
##   .. .. .. ..$ : chr [1:3] "y" "S" "D"
##   .. .. .. ..$ : chr [1:3] "S" "D" "S:D"
##   .. ..- attr(*, "term.labels")= chr [1:3] "S" "D" "S:D"
##   .. ..- attr(*, "order")= int [1:3] 1 1 2
##   .. ..- attr(*, "intercept")= int 1
##   .. ..- attr(*, "response")= int 1
##   .. ..- attr(*, ".Environment")=<environment: R_GlobalEnv> 
##   .. ..- attr(*, "predvars")= language list(y, S, D)
##   .. ..- attr(*, "dataClasses")= Named chr [1:3] "numeric" "character" "character"
##   .. .. ..- attr(*, "names")= chr [1:3] "y" "S" "D"
##  $ model        :'data.frame':   61 obs. of  3 variables:
##   ..$ y: int [1:61] 102 136 84 101 76 125 134 133 117 124 ...
##   ..$ S: chr [1:61] "WL" "WL" "WL" "WL" ...
##   ..$ D: chr [1:61] "RecipientPre" "RecipientPre" "RecipientPre" "RecipientPre" ...
##   ..- attr(*, "terms")=Classes 'terms', 'formula'  language y ~ S * D
##   .. .. ..- attr(*, "variables")= language list(y, S, D)
##   .. .. ..- attr(*, "factors")= int [1:3, 1:3] 0 1 0 0 0 1 0 1 1
##   .. .. .. ..- attr(*, "dimnames")=List of 2
##   .. .. .. .. ..$ : chr [1:3] "y" "S" "D"
##   .. .. .. .. ..$ : chr [1:3] "S" "D" "S:D"
##   .. .. ..- attr(*, "term.labels")= chr [1:3] "S" "D" "S:D"
##   .. .. ..- attr(*, "order")= int [1:3] 1 1 2
##   .. .. ..- attr(*, "intercept")= int 1
##   .. .. ..- attr(*, "response")= int 1
##   .. .. ..- attr(*, ".Environment")=<environment: R_GlobalEnv> 
##   .. .. ..- attr(*, "predvars")= language list(y, S, D)
##   .. .. ..- attr(*, "dataClasses")= Named chr [1:3] "numeric" "character" "character"
##   .. .. .. ..- attr(*, "names")= chr [1:3] "y" "S" "D"
##  - attr(*, "class")= chr "lm"
summary(DSRhizRich)
## 
## Call:
## lm(formula = y ~ S * D)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -42.533  -7.533   2.000   7.071  35.167 
## 
## Coefficients: (5 not defined because of singularities)
##                    Estimate Std. Error t value Pr(>|t|)    
## (Intercept)         125.400      6.941  18.066  < 2e-16 ***
## SHH                   8.600     12.986   0.662 0.510616    
## SPS                  -3.400     10.412  -0.327 0.745266    
## SWL                 -18.567      9.398  -1.976 0.053328 .  
## DRecipientPost       31.095      7.573   4.106 0.000137 ***
## DRecipientPre            NA         NA      NA       NA    
## SHH:DRecipientPost  -26.562     13.924  -1.908 0.061754 .  
## SPS:DRecipientPost  -18.962     11.560  -1.640 0.106764    
## SWL:DRecipientPost       NA         NA      NA       NA    
## SHH:DRecipientPre        NA         NA      NA       NA    
## SPS:DRecipientPre        NA         NA      NA       NA    
## SWL:DRecipientPre        NA         NA      NA       NA    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 15.52 on 54 degrees of freedom
## Multiple R-squared:  0.3019, Adjusted R-squared:  0.2243 
## F-statistic: 3.892 on 6 and 54 DF,  p-value: 0.002662
y <- ssuguildRICH$Rhizophilic_Richness
str(y)
##  int [1:61] 102 136 84 101 76 125 134 133 117 124 ...
D <- ssuguildRICH$Description
S <- ssuguildRICH$Site
length(S)
## [1] 61
length(D)
## [1] 61
`?`(lm)

glm(y ~ D)
## 
## Call:  glm(formula = y ~ D)
## 
## Coefficients:
##    (Intercept)  DRecipientPost   DRecipientPre  
##        125.400          11.441          -8.983  
## 
## Degrees of Freedom: 60 Total (i.e. Null);  58 Residual
## Null Deviance:       18630 
## Residual Deviance: 14470     AIC: 514.7
G_SRhizRich <- glm(y ~ S)
str(G_SRhizRich)
## List of 30
##  $ coefficients     : Named num [1:4] 125.4 12.6 6.18 3.2
##   ..- attr(*, "names")= chr [1:4] "(Intercept)" "SHH" "SPS" "SWL"
##  $ residuals        : Named num [1:61] -26.6 7.4 -44.6 -27.6 -52.6 ...
##   ..- attr(*, "names")= chr [1:61] "1" "2" "3" "4" ...
##  $ fitted.values    : Named num [1:61] 129 129 129 129 129 ...
##   ..- attr(*, "names")= chr [1:61] "1" "2" "3" "4" ...
##  $ effects          : Named num [1:61] -1030.1 29.7 11.9 -6.4 -44.7 ...
##   ..- attr(*, "names")= chr [1:61] "(Intercept)" "SHH" "SPS" "SWL" ...
##  $ R                : num [1:4, 1:4] -7.81 0 0 0 -2.18 ...
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : chr [1:4] "(Intercept)" "SHH" "SPS" "SWL"
##   .. ..$ : chr [1:4] "(Intercept)" "SHH" "SPS" "SWL"
##  $ rank             : int 4
##  $ qr               :List of 5
##   ..$ qr   : num [1:61, 1:4] -7.81 0.128 0.128 0.128 0.128 ...
##   .. ..- attr(*, "dimnames")=List of 2
##   .. .. ..$ : chr [1:61] "1" "2" "3" "4" ...
##   .. .. ..$ : chr [1:4] "(Intercept)" "SHH" "SPS" "SWL"
##   ..$ rank : int 4
##   ..$ qraux: num [1:4] 1.13 1.07 1.11 1.07
##   ..$ pivot: int [1:4] 1 2 3 4
##   ..$ tol  : num 1e-11
##   ..- attr(*, "class")= chr "qr"
##  $ family           :List of 11
##   ..$ family    : chr "gaussian"
##   ..$ link      : chr "identity"
##   ..$ linkfun   :function (mu)  
##   ..$ linkinv   :function (eta)  
##   ..$ variance  :function (mu)  
##   ..$ dev.resids:function (y, mu, wt)  
##   ..$ aic       :function (y, n, mu, wt, dev)  
##   ..$ mu.eta    :function (eta)  
##   ..$ initialize:  expression({  n <- rep.int(1, nobs)  if (is.null(etastart) && is.null(start) && is.null(mustart) &&  ((family$lin| __truncated__
##   ..$ validmu   :function (mu)  
##   ..$ valideta  :function (eta)  
##   ..- attr(*, "class")= chr "family"
##  $ linear.predictors: Named num [1:61] 129 129 129 129 129 ...
##   ..- attr(*, "names")= chr [1:61] "1" "2" "3" "4" ...
##  $ deviance         : num 17571
##  $ aic              : num 529
##  $ null.deviance    : num 18634
##  $ iter             : int 2
##  $ weights          : Named num [1:61] 1 1 1 1 1 1 1 1 1 1 ...
##   ..- attr(*, "names")= chr [1:61] "1" "2" "3" "4" ...
##  $ prior.weights    : Named num [1:61] 1 1 1 1 1 1 1 1 1 1 ...
##   ..- attr(*, "names")= chr [1:61] "1" "2" "3" "4" ...
##  $ df.residual      : int 57
##  $ df.null          : int 60
##  $ y                : Named int [1:61] 102 136 84 101 76 125 134 133 117 124 ...
##   ..- attr(*, "names")= chr [1:61] "1" "2" "3" "4" ...
##  $ converged        : logi TRUE
##  $ boundary         : logi FALSE
##  $ model            :'data.frame':   61 obs. of  2 variables:
##   ..$ y: int [1:61] 102 136 84 101 76 125 134 133 117 124 ...
##   ..$ S: chr [1:61] "WL" "WL" "WL" "WL" ...
##   ..- attr(*, "terms")=Classes 'terms', 'formula'  language y ~ S
##   .. .. ..- attr(*, "variables")= language list(y, S)
##   .. .. ..- attr(*, "factors")= int [1:2, 1] 0 1
##   .. .. .. ..- attr(*, "dimnames")=List of 2
##   .. .. .. .. ..$ : chr [1:2] "y" "S"
##   .. .. .. .. ..$ : chr "S"
##   .. .. ..- attr(*, "term.labels")= chr "S"
##   .. .. ..- attr(*, "order")= int 1
##   .. .. ..- attr(*, "intercept")= int 1
##   .. .. ..- attr(*, "response")= int 1
##   .. .. ..- attr(*, ".Environment")=<environment: R_GlobalEnv> 
##   .. .. ..- attr(*, "predvars")= language list(y, S)
##   .. .. ..- attr(*, "dataClasses")= Named chr [1:2] "numeric" "character"
##   .. .. .. ..- attr(*, "names")= chr [1:2] "y" "S"
##  $ call             : language glm(formula = y ~ S)
##  $ formula          :Class 'formula'  language y ~ S
##   .. ..- attr(*, ".Environment")=<environment: R_GlobalEnv> 
##  $ terms            :Classes 'terms', 'formula'  language y ~ S
##   .. ..- attr(*, "variables")= language list(y, S)
##   .. ..- attr(*, "factors")= int [1:2, 1] 0 1
##   .. .. ..- attr(*, "dimnames")=List of 2
##   .. .. .. ..$ : chr [1:2] "y" "S"
##   .. .. .. ..$ : chr "S"
##   .. ..- attr(*, "term.labels")= chr "S"
##   .. ..- attr(*, "order")= int 1
##   .. ..- attr(*, "intercept")= int 1
##   .. ..- attr(*, "response")= int 1
##   .. ..- attr(*, ".Environment")=<environment: R_GlobalEnv> 
##   .. ..- attr(*, "predvars")= language list(y, S)
##   .. ..- attr(*, "dataClasses")= Named chr [1:2] "numeric" "character"
##   .. .. ..- attr(*, "names")= chr [1:2] "y" "S"
##  $ data             :<environment: R_GlobalEnv> 
##  $ offset           : NULL
##  $ control          :List of 3
##   ..$ epsilon: num 1e-08
##   ..$ maxit  : num 25
##   ..$ trace  : logi FALSE
##  $ method           : chr "glm.fit"
##  $ contrasts        :List of 1
##   ..$ S: chr "contr.treatment"
##  $ xlevels          :List of 1
##   ..$ S: chr [1:4] "DS" "HH" "PS" "WL"
##  - attr(*, "class")= chr [1:2] "glm" "lm"
G_DRhizRich <- glm(y ~ D)
str(G_DRhizRich)
## List of 30
##  $ coefficients     : Named num [1:3] 125.4 11.44 -8.98
##   ..- attr(*, "names")= chr [1:3] "(Intercept)" "DRecipientPost" "DRecipientPre"
##  $ residuals        : Named num [1:61] -14.4 19.6 -32.4 -15.4 -40.4 ...
##   ..- attr(*, "names")= chr [1:61] "1" "2" "3" "4" ...
##  $ fitted.values    : Named num [1:61] 116 116 116 116 116 ...
##   ..- attr(*, "names")= chr [1:61] "1" "2" "3" "4" ...
##  $ effects          : Named num [1:61] -1030.1 62.3 16.9 -13.5 -38.5 ...
##   ..- attr(*, "names")= chr [1:61] "(Intercept)" "DRecipientPost" "DRecipientPre" "" ...
##  $ R                : num [1:3, 1:3] -7.81 0 0 -5.63 3.5 ...
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : chr [1:3] "(Intercept)" "DRecipientPost" "DRecipientPre"
##   .. ..$ : chr [1:3] "(Intercept)" "DRecipientPost" "DRecipientPre"
##  $ rank             : int 3
##  $ qr               :List of 5
##   ..$ qr   : num [1:61, 1:3] -7.81 0.128 0.128 0.128 0.128 ...
##   .. ..- attr(*, "dimnames")=List of 2
##   .. .. ..$ : chr [1:61] "1" "2" "3" "4" ...
##   .. .. ..$ : chr [1:3] "(Intercept)" "DRecipientPost" "DRecipientPre"
##   ..$ rank : int 3
##   ..$ qraux: num [1:3] 1.13 1.18 1.12
##   ..$ pivot: int [1:3] 1 2 3
##   ..$ tol  : num 1e-11
##   ..- attr(*, "class")= chr "qr"
##  $ family           :List of 11
##   ..$ family    : chr "gaussian"
##   ..$ link      : chr "identity"
##   ..$ linkfun   :function (mu)  
##   ..$ linkinv   :function (eta)  
##   ..$ variance  :function (mu)  
##   ..$ dev.resids:function (y, mu, wt)  
##   ..$ aic       :function (y, n, mu, wt, dev)  
##   ..$ mu.eta    :function (eta)  
##   ..$ initialize:  expression({  n <- rep.int(1, nobs)  if (is.null(etastart) && is.null(start) && is.null(mustart) &&  ((family$lin| __truncated__
##   ..$ validmu   :function (mu)  
##   ..$ valideta  :function (eta)  
##   ..- attr(*, "class")= chr "family"
##  $ linear.predictors: Named num [1:61] 116 116 116 116 116 ...
##   ..- attr(*, "names")= chr [1:61] "1" "2" "3" "4" ...
##  $ deviance         : num 14472
##  $ aic              : num 515
##  $ null.deviance    : num 18634
##  $ iter             : int 2
##  $ weights          : Named num [1:61] 1 1 1 1 1 1 1 1 1 1 ...
##   ..- attr(*, "names")= chr [1:61] "1" "2" "3" "4" ...
##  $ prior.weights    : Named num [1:61] 1 1 1 1 1 1 1 1 1 1 ...
##   ..- attr(*, "names")= chr [1:61] "1" "2" "3" "4" ...
##  $ df.residual      : int 58
##  $ df.null          : int 60
##  $ y                : Named int [1:61] 102 136 84 101 76 125 134 133 117 124 ...
##   ..- attr(*, "names")= chr [1:61] "1" "2" "3" "4" ...
##  $ converged        : logi TRUE
##  $ boundary         : logi FALSE
##  $ model            :'data.frame':   61 obs. of  2 variables:
##   ..$ y: int [1:61] 102 136 84 101 76 125 134 133 117 124 ...
##   ..$ D: chr [1:61] "RecipientPre" "RecipientPre" "RecipientPre" "RecipientPre" ...
##   ..- attr(*, "terms")=Classes 'terms', 'formula'  language y ~ D
##   .. .. ..- attr(*, "variables")= language list(y, D)
##   .. .. ..- attr(*, "factors")= int [1:2, 1] 0 1
##   .. .. .. ..- attr(*, "dimnames")=List of 2
##   .. .. .. .. ..$ : chr [1:2] "y" "D"
##   .. .. .. .. ..$ : chr "D"
##   .. .. ..- attr(*, "term.labels")= chr "D"
##   .. .. ..- attr(*, "order")= int 1
##   .. .. ..- attr(*, "intercept")= int 1
##   .. .. ..- attr(*, "response")= int 1
##   .. .. ..- attr(*, ".Environment")=<environment: R_GlobalEnv> 
##   .. .. ..- attr(*, "predvars")= language list(y, D)
##   .. .. ..- attr(*, "dataClasses")= Named chr [1:2] "numeric" "character"
##   .. .. .. ..- attr(*, "names")= chr [1:2] "y" "D"
##  $ call             : language glm(formula = y ~ D)
##  $ formula          :Class 'formula'  language y ~ D
##   .. ..- attr(*, ".Environment")=<environment: R_GlobalEnv> 
##  $ terms            :Classes 'terms', 'formula'  language y ~ D
##   .. ..- attr(*, "variables")= language list(y, D)
##   .. ..- attr(*, "factors")= int [1:2, 1] 0 1
##   .. .. ..- attr(*, "dimnames")=List of 2
##   .. .. .. ..$ : chr [1:2] "y" "D"
##   .. .. .. ..$ : chr "D"
##   .. ..- attr(*, "term.labels")= chr "D"
##   .. ..- attr(*, "order")= int 1
##   .. ..- attr(*, "intercept")= int 1
##   .. ..- attr(*, "response")= int 1
##   .. ..- attr(*, ".Environment")=<environment: R_GlobalEnv> 
##   .. ..- attr(*, "predvars")= language list(y, D)
##   .. ..- attr(*, "dataClasses")= Named chr [1:2] "numeric" "character"
##   .. .. ..- attr(*, "names")= chr [1:2] "y" "D"
##  $ data             :<environment: R_GlobalEnv> 
##  $ offset           : NULL
##  $ control          :List of 3
##   ..$ epsilon: num 1e-08
##   ..$ maxit  : num 25
##   ..$ trace  : logi FALSE
##  $ method           : chr "glm.fit"
##  $ contrasts        :List of 1
##   ..$ D: chr "contr.treatment"
##  $ xlevels          :List of 1
##   ..$ D: chr [1:3] "Donor" "RecipientPost" "RecipientPre"
##  - attr(*, "class")= chr [1:2] "glm" "lm"
summary(G_DRhizRich)
## 
## Call:
## glm(formula = y ~ D)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -40.841   -8.400    3.159   10.159   26.159  
## 
## Coefficients:
##                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)     125.400      7.064  17.751   <2e-16 ***
## DRecipientPost   11.441      7.455   1.535     0.13    
## DRecipientPre    -8.983      8.408  -1.068     0.29    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for gaussian family taken to be 249.5173)
## 
##     Null deviance: 18634  on 60  degrees of freedom
## Residual deviance: 14472  on 58  degrees of freedom
## AIC: 514.73
## 
## Number of Fisher Scoring iterations: 2
DSRhizRich <- lm(y ~ S * D)
str(DSRhizRich)
## List of 13
##  $ coefficients : Named num [1:12] 125.4 8.6 -3.4 -18.6 31.1 ...
##   ..- attr(*, "names")= chr [1:12] "(Intercept)" "SHH" "SPS" "SWL" ...
##  $ residuals    : Named num [1:61] -4.83 29.17 -22.83 -5.83 -30.83 ...
##   ..- attr(*, "names")= chr [1:61] "1" "2" "3" "4" ...
##  $ effects      : Named num [1:61] -1030.1 29.7 11.9 -6.4 58.6 ...
##   ..- attr(*, "names")= chr [1:61] "(Intercept)" "SHH" "SPS" "SWL" ...
##  $ rank         : int 7
##  $ fitted.values: Named num [1:61] 107 107 107 107 107 ...
##   ..- attr(*, "names")= chr [1:61] "1" "2" "3" "4" ...
##  $ assign       : int [1:12] 0 1 1 1 2 2 3 3 3 3 ...
##  $ qr           :List of 5
##   ..$ qr   : num [1:61, 1:12] -7.81 0.128 0.128 0.128 0.128 ...
##   .. ..- attr(*, "dimnames")=List of 2
##   .. .. ..$ : chr [1:61] "1" "2" "3" "4" ...
##   .. .. ..$ : chr [1:12] "(Intercept)" "SHH" "SPS" "SWL" ...
##   .. ..- attr(*, "assign")= int [1:12] 0 1 1 1 2 2 3 3 3 3 ...
##   .. ..- attr(*, "contrasts")=List of 2
##   .. .. ..$ S: chr "contr.treatment"
##   .. .. ..$ D: chr "contr.treatment"
##   ..$ qraux: num [1:12] 1.13 1.07 1.11 1.07 1.16 ...
##   ..$ pivot: int [1:12] 1 2 3 4 5 7 8 6 9 10 ...
##   ..$ tol  : num 1e-07
##   ..$ rank : int 7
##   ..- attr(*, "class")= chr "qr"
##  $ df.residual  : int 54
##  $ contrasts    :List of 2
##   ..$ S: chr "contr.treatment"
##   ..$ D: chr "contr.treatment"
##  $ xlevels      :List of 2
##   ..$ S: chr [1:4] "DS" "HH" "PS" "WL"
##   ..$ D: chr [1:3] "Donor" "RecipientPost" "RecipientPre"
##  $ call         : language lm(formula = y ~ S * D)
##  $ terms        :Classes 'terms', 'formula'  language y ~ S * D
##   .. ..- attr(*, "variables")= language list(y, S, D)
##   .. ..- attr(*, "factors")= int [1:3, 1:3] 0 1 0 0 0 1 0 1 1
##   .. .. ..- attr(*, "dimnames")=List of 2
##   .. .. .. ..$ : chr [1:3] "y" "S" "D"
##   .. .. .. ..$ : chr [1:3] "S" "D" "S:D"
##   .. ..- attr(*, "term.labels")= chr [1:3] "S" "D" "S:D"
##   .. ..- attr(*, "order")= int [1:3] 1 1 2
##   .. ..- attr(*, "intercept")= int 1
##   .. ..- attr(*, "response")= int 1
##   .. ..- attr(*, ".Environment")=<environment: R_GlobalEnv> 
##   .. ..- attr(*, "predvars")= language list(y, S, D)
##   .. ..- attr(*, "dataClasses")= Named chr [1:3] "numeric" "character" "character"
##   .. .. ..- attr(*, "names")= chr [1:3] "y" "S" "D"
##  $ model        :'data.frame':   61 obs. of  3 variables:
##   ..$ y: int [1:61] 102 136 84 101 76 125 134 133 117 124 ...
##   ..$ S: chr [1:61] "WL" "WL" "WL" "WL" ...
##   ..$ D: chr [1:61] "RecipientPre" "RecipientPre" "RecipientPre" "RecipientPre" ...
##   ..- attr(*, "terms")=Classes 'terms', 'formula'  language y ~ S * D
##   .. .. ..- attr(*, "variables")= language list(y, S, D)
##   .. .. ..- attr(*, "factors")= int [1:3, 1:3] 0 1 0 0 0 1 0 1 1
##   .. .. .. ..- attr(*, "dimnames")=List of 2
##   .. .. .. .. ..$ : chr [1:3] "y" "S" "D"
##   .. .. .. .. ..$ : chr [1:3] "S" "D" "S:D"
##   .. .. ..- attr(*, "term.labels")= chr [1:3] "S" "D" "S:D"
##   .. .. ..- attr(*, "order")= int [1:3] 1 1 2
##   .. .. ..- attr(*, "intercept")= int 1
##   .. .. ..- attr(*, "response")= int 1
##   .. .. ..- attr(*, ".Environment")=<environment: R_GlobalEnv> 
##   .. .. ..- attr(*, "predvars")= language list(y, S, D)
##   .. .. ..- attr(*, "dataClasses")= Named chr [1:3] "numeric" "character" "character"
##   .. .. .. ..- attr(*, "names")= chr [1:3] "y" "S" "D"
##  - attr(*, "class")= chr "lm"
summary(DSRhizRich)
## 
## Call:
## lm(formula = y ~ S * D)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -42.533  -7.533   2.000   7.071  35.167 
## 
## Coefficients: (5 not defined because of singularities)
##                    Estimate Std. Error t value Pr(>|t|)    
## (Intercept)         125.400      6.941  18.066  < 2e-16 ***
## SHH                   8.600     12.986   0.662 0.510616    
## SPS                  -3.400     10.412  -0.327 0.745266    
## SWL                 -18.567      9.398  -1.976 0.053328 .  
## DRecipientPost       31.095      7.573   4.106 0.000137 ***
## DRecipientPre            NA         NA      NA       NA    
## SHH:DRecipientPost  -26.562     13.924  -1.908 0.061754 .  
## SPS:DRecipientPost  -18.962     11.560  -1.640 0.106764    
## SWL:DRecipientPost       NA         NA      NA       NA    
## SHH:DRecipientPre        NA         NA      NA       NA    
## SPS:DRecipientPre        NA         NA      NA       NA    
## SWL:DRecipientPre        NA         NA      NA       NA    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 15.52 on 54 degrees of freedom
## Multiple R-squared:  0.3019, Adjusted R-squared:  0.2243 
## F-statistic: 3.892 on 6 and 54 DF,  p-value: 0.002662
lm(y ~ D * S)
## 
## Call:
## lm(formula = y ~ D * S)
## 
## Coefficients:
##        (Intercept)      DRecipientPost       DRecipientPre  
##             125.40               12.53              -18.57  
##                SHH                 SPS                 SWL  
##              27.17               15.17                  NA  
## DRecipientPost:SHH   DRecipientPre:SHH  DRecipientPost:SPS  
##             -26.56                  NA              -18.96  
##  DRecipientPre:SPS  DRecipientPost:SWL   DRecipientPre:SWL  
##                 NA                  NA                  NA
lm(y ~ D + S)
## 
## Call:
## lm(formula = y ~ D + S)
## 
## Coefficients:
##    (Intercept)  DRecipientPost   DRecipientPre             SHH  
##        125.400           9.018         -10.375           5.864  
##            SPS             SWL  
##          1.244              NA
`?`(glm)
glm(y ~ D)
## 
## Call:  glm(formula = y ~ D)
## 
## Coefficients:
##    (Intercept)  DRecipientPost   DRecipientPre  
##        125.400          11.441          -8.983  
## 
## Degrees of Freedom: 60 Total (i.e. Null);  58 Residual
## Null Deviance:       18630 
## Residual Deviance: 14470     AIC: 514.7
# Description+Site

start ‘2017_16_07_SSU_Boxplots.R’ Plotting

# richness
plot.df <- ssuguildREADRICH
str(plot.df)
## 'data.frame':    183 obs. of  16 variables:
##  $ ID                   : chr  "MM933M" "MM933M" "MM933M" "MM936M" ...
##  $ Functional.Group     : chr  "Ancestral" "Edaphophilic" "Rhizophilic" "Ancestral" ...
##  $ OTU_Richness_Sample  : int  26 2 102 20 1 84 19 1 101 20 ...
##  $ Read_Abundance_Sample: num  194 23.2 616.6 103.9 12.2 ...
##  $ #SampleID            : int  5 5 5 3 3 3 1 1 1 2 ...
##  $ BarcodeSequence      : chr  "ATCCCGTATCGATTGG" "ATCCCGTATCGATTGG" "ATCCCGTATCGATTGG" "GCAACCTTAGAGTGTG" ...
##  $ LinkerPrimerSequence : logi  NA NA NA NA NA NA ...
##  $ Location             : chr  "E1" "E1" "E1" "C1" ...
##  $ Project              : chr  "MM-Salvage" "MM-Salvage" "MM-Salvage" "MM-Salvage" ...
##  $ Year                 : chr  "2015" "2015" "2015" "2015" ...
##  $ Site                 : chr  "WL" "WL" "WL" "WL" ...
##  $ Treat                : Factor w/ 6 levels "C","D","F","O",..: 4 4 4 4 4 4 4 4 4 4 ...
##  $ Rep                  : Factor w/ 6 levels "1","2","3","4",..: 5 5 5 3 3 3 1 1 1 2 ...
##  $ HL                   : num  0.264 0.264 0.264 0.264 0.264 ...
##  $ TSDel                : chr  NA NA NA NA ...
##  $ Description          : chr  "RecipientPre" "RecipientPre" "RecipientPre" "RecipientPre" ...
richbox <- ggplot(plot.df, aes(x = Functional.Group, y = OTU_Richness_Sample)) + 
    geom_boxplot(aes(fill = Description))
richbox

richbox <- richbox + theme_bw(base_size = 15) + xlab("Functional Group") + ylab("AMF Taxa Richness")
richbox

richbox <- richbox + scale_fill_manual(values = c("red", "darkslateblue", "gold3"))
richbox

## reads
plot.read <- ssuguildREADRICH
str(plot.read)

reachbox <- ggplot(plot.read, aes(x = Functional.Group, y = OTU_Richness_Sample)) + 
    geom_boxplot(aes(fill = Description))
reachbox

reachbox <- reachbox + theme_bw(base_size = 15) + xlab("Functional Group") + 
    ylab("AMF Taxa Reads")
reachbox

reachbox <- reachbox + scale_fill_manual(values = c("red", "darkslateblue", 
    "gold3"))
reachbox

ADAPTED FROM #GLMM for roots and soil OTU Richness SSU by functional group-modified from “GLM_NAU_6_2017” #—different models for each functional group #Author: Michala Phillips #MODIFY FOR THIS DATA

library(car)
require(MASS)
library(lme4)

Choose distribution based on data

find the right distribution for data

ssuguildREADRICHlong$OTU_Richness_Sample.t <- ssuguildREADRICHlong$OTU_Richness_Sample + 
    1
qqp(ssuguildREADRICHlong$OTU_Richness_Sample.t, "norm")

## [1] 5 3
# LOGNORMAL
qqp(ssuguildREADRICHlong$OTU_Richness_Sample, "lnorm")

## [1] 5 3
# NEG BINOMIAL
nbinom <- fitdistr(ssuguildREADRICHlong$OTU_Richness_Sample.t, "Negative Binomial")
qqp(ssuguildREADRICHlong$OTU_Richness_Sample.t, "nbinom", size = nbinom$estimate[[1]], 
    mu = nbinom$estimate[[2]])

## [1] 5 3
# POISSON
poisson <- fitdistr(ssuguildREADRICHlong$OTU_Richness_Sample.t, "Poisson")
qqp(ssuguildREADRICHlong$OTU_Richness_Sample.t, "pois", lambda = poisson$estimate)

## [1] 5 3
# GAMMA
gamma <- fitdistr(ssuguildREADRICHlong$OTU_Richness_Sample.t, "gamma")
qqp(ssuguildREADRICHlong$OTU_Richness_Sample.t, "gamma", shape = gamma$estimate[[1]], 
    rate = gamma$estimate[[2]])

## [1] 5 3
# 
modelingAMF <- glm(OTU_Richness_Sample.t ~ Site + Treat, data = ssuguildREADRICHlong, 
    family = gaussian(link = "identity"))
modelingAMF
## 
## Call:  glm(formula = OTU_Richness_Sample.t ~ Site + Treat, family = gaussian(link = "identity"), 
##     data = ssuguildREADRICHlong)
## 
## Coefficients:
## (Intercept)       SiteHH       SitePS       SiteWL       TreatD  
##    174.4580      -3.4348     -10.8522      -8.4202     -12.8889  
##      TreatF       TreatO       TreatS       TreatT  
##     11.1111     -24.0580      -0.8703       1.2222  
## 
## Degrees of Freedom: 60 Total (i.e. Null);  52 Residual
## Null Deviance:       26340 
## Residual Deviance: 16940     AIC: 536.3
modelingAMF_MP <- glm(OTU_Richness_Sample.t ~ Site + Treat + Description + Rep + 
    Year, data = ssuguildREADRICHlong, family = gaussian(link = "identity"))
plot(modelingAMF_MP)

stepAIC(modelingAMF_MP)
## Start:  AIC=519.78
## OTU_Richness_Sample.t ~ Site + Treat + Description + Rep + Year
## 
## 
## Step:  AIC=519.78
## OTU_Richness_Sample.t ~ Site + Treat + Description + Rep
## 
## 
## Step:  AIC=519.78
## OTU_Richness_Sample.t ~ Site + Treat + Rep
## 
##         Df Deviance    AIC
## - Site   3    12025 519.43
## <none>        10961 519.78
## - Rep    5    16942 536.34
## - Treat  5    20352 547.52
## 
## Step:  AIC=519.43
## OTU_Richness_Sample.t ~ Treat + Rep
## 
##         Df Deviance    AIC
## <none>        12025 519.43
## - Rep    5    17693 532.98
## - Treat  5    21477 544.81
## 
## Call:  glm(formula = OTU_Richness_Sample.t ~ Treat + Rep, family = gaussian(link = "identity"), 
##     data = ssuguildREADRICHlong)
## 
## Coefficients:
## (Intercept)       TreatD       TreatF       TreatO       TreatS  
##     165.889      -12.889       11.111      -28.319       -1.981  
##      TreatT         Rep2         Rep3         Rep4         Rep5  
##       1.222       11.737       -8.737       21.431       -1.569  
##        Rep6  
##      35.431  
## 
## Degrees of Freedom: 60 Total (i.e. Null);  50 Residual
## Null Deviance:       26340 
## Residual Deviance: 12020     AIC: 519.4
modelingAMF_MP
## 
## Call:  glm(formula = OTU_Richness_Sample.t ~ Site + Treat + Description + 
##     Rep + Year, family = gaussian(link = "identity"), data = ssuguildREADRICHlong)
## 
## Coefficients:
##              (Intercept)                    SiteHH  
##                  177.106                    -6.603  
##                   SitePS                    SiteWL  
##                  -13.930                   -13.110  
##                   TreatD                    TreatF  
##                  -12.889                    11.111  
##                   TreatO                    TreatS  
##                  -30.959                    -2.219  
##                   TreatT  DescriptionRecipientPost  
##                    1.222                        NA  
##  DescriptionRecipientPre                      Rep2  
##                       NA                    11.737  
##                     Rep3                      Rep4  
##                   -8.746                    21.866  
##                     Rep5                      Rep6  
##                   -3.592                    39.963  
##                 Year2017  
##                       NA  
## 
## Degrees of Freedom: 60 Total (i.e. Null);  47 Residual
## Null Deviance:       26340 
## Residual Deviance: 10960     AIC: 519.8
stepAIC(modelingAMF)
## Start:  AIC=536.34
## OTU_Richness_Sample.t ~ Site + Treat
## 
##         Df Deviance    AIC
## - Site   3    17693 532.98
## <none>        16942 536.34
## - Treat  5    24729 549.41
## 
## Step:  AIC=532.98
## OTU_Richness_Sample.t ~ Treat
## 
##         Df Deviance    AIC
## <none>        17693 532.98
## - Treat  5    26339 547.25
## 
## Call:  glm(formula = OTU_Richness_Sample.t ~ Treat, family = gaussian(link = "identity"), 
##     data = ssuguildREADRICHlong)
## 
## Coefficients:
## (Intercept)       TreatD       TreatF       TreatO       TreatS  
##    166.8889     -12.8889      11.1111     -22.4183      -0.7639  
##      TreatT  
##      1.2222  
## 
## Degrees of Freedom: 60 Total (i.e. Null);  55 Residual
## Null Deviance:       26340 
## Residual Deviance: 17690     AIC: 533
modelingAMFfinal <- glm(OTU_Richness_Sample.t ~ Site + Treat, family = gaussian(link = "identity"), 
    data = ssuguildREADRICHlong)

summary(modelingAMFfinal)
## 
## Call:
## glm(formula = OTU_Richness_Sample.t ~ Site + Treat, family = gaussian(link = "identity"), 
##     data = ssuguildREADRICHlong)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -41.980   -8.400    3.265    9.740   32.962  
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 174.4580    11.3814  15.328  < 2e-16 ***
## SiteHH       -3.4348    10.5906  -0.324  0.74699    
## SitePS      -10.8522    10.2271  -1.061  0.29354    
## SiteWL       -8.4202     9.9487  -0.846  0.40123    
## TreatD      -12.8889     8.5090  -1.515  0.13590    
## TreatF       11.1111     8.5090   1.306  0.19737    
## TreatO      -24.0580     8.0233  -2.999  0.00415 ** 
## TreatS       -0.8703     8.7812  -0.099  0.92143    
## TreatT        1.2222     8.5090   0.144  0.88634    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for gaussian family taken to be 325.8149)
## 
##     Null deviance: 26339  on 60  degrees of freedom
## Residual deviance: 16942  on 52  degrees of freedom
## AIC: 536.34
## 
## Number of Fisher Scoring iterations: 2
plot(modelingAMFfinal)

modelingAMF <- glm(OTU_Richness_Sample.t ~ Site + Treat, data = ssuguildREADRICHlong, 
    family = gaussian(link = "identity"))
modelingAMF
## 
## Call:  glm(formula = OTU_Richness_Sample.t ~ Site + Treat, family = gaussian(link = "identity"), 
##     data = ssuguildREADRICHlong)
## 
## Coefficients:
## (Intercept)       SiteHH       SitePS       SiteWL       TreatD  
##    174.4580      -3.4348     -10.8522      -8.4202     -12.8889  
##      TreatF       TreatO       TreatS       TreatT  
##     11.1111     -24.0580      -0.8703       1.2222  
## 
## Degrees of Freedom: 60 Total (i.e. Null);  52 Residual
## Null Deviance:       26340 
## Residual Deviance: 16940     AIC: 536.3
stepAIC(modelingAMF)
## Start:  AIC=536.34
## OTU_Richness_Sample.t ~ Site + Treat
## 
##         Df Deviance    AIC
## - Site   3    17693 532.98
## <none>        16942 536.34
## - Treat  5    24729 549.41
## 
## Step:  AIC=532.98
## OTU_Richness_Sample.t ~ Treat
## 
##         Df Deviance    AIC
## <none>        17693 532.98
## - Treat  5    26339 547.25
## 
## Call:  glm(formula = OTU_Richness_Sample.t ~ Treat, family = gaussian(link = "identity"), 
##     data = ssuguildREADRICHlong)
## 
## Coefficients:
## (Intercept)       TreatD       TreatF       TreatO       TreatS  
##    166.8889     -12.8889      11.1111     -22.4183      -0.7639  
##      TreatT  
##      1.2222  
## 
## Degrees of Freedom: 60 Total (i.e. Null);  55 Residual
## Null Deviance:       26340 
## Residual Deviance: 17690     AIC: 533
modelingAMFfinal <- glm(OTU_Richness_Sample.t ~ Site + Treat, family = gaussian(link = "identity"), 
    data = ssuguildREADRICHlong)

summary(modelingAMFfinal)
## 
## Call:
## glm(formula = OTU_Richness_Sample.t ~ Site + Treat, family = gaussian(link = "identity"), 
##     data = ssuguildREADRICHlong)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -41.980   -8.400    3.265    9.740   32.962  
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 174.4580    11.3814  15.328  < 2e-16 ***
## SiteHH       -3.4348    10.5906  -0.324  0.74699    
## SitePS      -10.8522    10.2271  -1.061  0.29354    
## SiteWL       -8.4202     9.9487  -0.846  0.40123    
## TreatD      -12.8889     8.5090  -1.515  0.13590    
## TreatF       11.1111     8.5090   1.306  0.19737    
## TreatO      -24.0580     8.0233  -2.999  0.00415 ** 
## TreatS       -0.8703     8.7812  -0.099  0.92143    
## TreatT        1.2222     8.5090   0.144  0.88634    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for gaussian family taken to be 325.8149)
## 
##     Null deviance: 26339  on 60  degrees of freedom
## Residual deviance: 16942  on 52  degrees of freedom
## AIC: 536.34
## 
## Number of Fisher Scoring iterations: 2
plot(modelingAMFfinal)

#Working on examining other variables

modelingAMF <- glm(OTU_Richness_Sample.t ~ Treat + Site, data = ssuguildREADRICHlong, 
    family = gaussian(link = "identity"))

Plotting

plot.df <- ssuguildREADRICH
str(plot.df)
## 'data.frame':    183 obs. of  16 variables:
##  $ ID                   : chr  "MM933M" "MM933M" "MM933M" "MM936M" ...
##  $ Functional.Group     : chr  "Ancestral" "Edaphophilic" "Rhizophilic" "Ancestral" ...
##  $ OTU_Richness_Sample  : int  26 2 102 20 1 84 19 1 101 20 ...
##  $ Read_Abundance_Sample: num  194 23.2 616.6 103.9 12.2 ...
##  $ #SampleID            : int  5 5 5 3 3 3 1 1 1 2 ...
##  $ BarcodeSequence      : chr  "ATCCCGTATCGATTGG" "ATCCCGTATCGATTGG" "ATCCCGTATCGATTGG" "GCAACCTTAGAGTGTG" ...
##  $ LinkerPrimerSequence : logi  NA NA NA NA NA NA ...
##  $ Location             : chr  "E1" "E1" "E1" "C1" ...
##  $ Project              : chr  "MM-Salvage" "MM-Salvage" "MM-Salvage" "MM-Salvage" ...
##  $ Year                 : chr  "2015" "2015" "2015" "2015" ...
##  $ Site                 : chr  "WL" "WL" "WL" "WL" ...
##  $ Treat                : Factor w/ 6 levels "C","D","F","O",..: 4 4 4 4 4 4 4 4 4 4 ...
##  $ Rep                  : Factor w/ 6 levels "1","2","3","4",..: 5 5 5 3 3 3 1 1 1 2 ...
##  $ HL                   : num  0.264 0.264 0.264 0.264 0.264 ...
##  $ TSDel                : chr  NA NA NA NA ...
##  $ Description          : chr  "RecipientPre" "RecipientPre" "RecipientPre" "RecipientPre" ...
richbox <- ggplot(plot.df, aes(x = Functional.Group, y = OTU_Richness_Sample)) + 
    geom_boxplot(aes(fill = Description))
richbox

richbox <- richbox + theme_bw(base_size = 15) + xlab("Functional Group") + ylab("AMF Taxa Richness")
richbox

richbox <- richbox + scale_fill_manual(values = c("red", "darkslateblue", "gold3"))
richbox

ADAPTED FROM #GLMM for roots and soil OTU Richness SSU by functional group-modified from “GLM_NAU_6_2017” #—different models for each functional group #Author: Michala Phillips #MODIFY FOR THIS DATA #Tried with negative binomial but it didn’t work

# salvage.div<-adonis2(bird.dist~DIVERSITY, data=birds, permutations = 999,
# method='bray', strata='PLOT') bird.div